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SUMMARY 


Reiner’ has observed a high pressure amounting to half an 
atmosphere at the center of the airspace between two parallel 
plates when they are very close together and one of them rotates 
These observations were interpreted by him as being due to as 
sumed non-Newtonian properties in air. Since aerodvnamic 
theory is founded on the Navier-Stokes equations which do not 
allow for non-Newtonian effects, it seemed desirable to explore 
possible conditions in which these equations might predict a high 
pressure between the plates. The effect of geometrical and dy- 
namical imperfections in the experimental apparatus is calcu 


lated, two particular cases being considered in detail \) asmall 


error in the perpendicularity of either of the plates to the axis of 


rotation and (B) a small vibration of the rotor parallel to its axis 
When the fluid is assumed incompressible, the mean value of th« 
excess pressure is found to be zero. The compressibility of the 
air changes this result, and mean pressures of the same order of 
magnitude as those observed experimentally are then predicted 
rhe effect is compared with that of compressibility in ait 
lubricated bearings \ corollary of the analysis leads to the 
conclusion that, when two wavy or rough surfaces slide over onc 
another, the effect of compressibility is to make the air laver be 


tween them exert a pressure tending to keep them apart 


INTRODUCTION 


ne AND REINER’ describe a most interesting 
experiment which was demonstrated by Professor 
Reiner at the 9th International Congress of Applied 
Mechanics at Brussels, 1956. A flat rotor dise 6.7 cm. 


in diameter was spun at about 7,000 r.p.m. about its 
axis of symmetry opposite a fixed stator disc which was 
also perpendicular to this axis. The stator was then 
moved slowly toward the rotor by means of a fine 
micrometer screw. A hole about 4 mm. in diameter 
in the center of the stator was connected to one side of a 
manometer, the other side of which was open to the 
atmosphere. When the gap between the plates was 
0.5 mm. apart, the manometer recorded a suction of 
about | cm. of water. This is close to what would be 


Presented at the Aerodynamics—I Session, Twenty-Fifth An 
nual Meeting, IAS, New York, January 28-31, 1957. Revised 


and received April 14, 1957 


expected, for Stewartson® has calculated that the suc 
tion at the center should be (3 (20) pw*a’, provided that 
the Reynolds Number pw? uw is not too large com 
pared with 1, where a is the radius of the discs, / the 
gap between them, w the angular velocity, and p and u 
the density and viscosity of air. This is 1.1 em. of 
water under the conditions of the above experiment 

This gap was then carefully closed, and the suction 
grew less till at 0.02 mm. —1.e., 20 microns—it disap 
peared. As the gap closed still further, the suction 
turned into a positive excess pressure. Excess pressures 
of up to half an atmosphere were reported before con 
tact of the two plates occurred. 

Reiner! interprets this experimental result as being 
due to a non-Newtonian or viscoelastic property of air 
analogous to that exhibited by certain fluids like rubber 
solution. These fluids have an elastic resistance to 
shear as well as viscosity and therefore possess a time 
of relaxation which manifests itself as a decreasing 
stress when the moving fluid is brought to rest or by 
the presence of ‘“‘cross stresses’ i.e., normal stresses 
in the direction of the velocity and velocity gradient 
when the fluid is sheared. Weissenberg* has demon 
strated the existence of such non-Newtonian properties 
with apparatus very similar to that described above 
When the space between rotor and stator was filled 
with certain fluids possessing long molecular chains, 
the pressure at the center was found to be higher than 
at the circumference. 

Reiner puts forward his experiment as a demonstra 
tion that air possesses elasticity for shear stress. How 
ever, the kinetic theory of gases predicts a time of re- 
laxation, u/po, for air equal to 1.9 K 107 sec., where 
po is atmospheric pressure; this time is, in fact, equal 
to the mean time between molecular collisions. Reiner 
cites Jeffreys as pointing out that it seems unlikely that 
any experiment could be devised which would involve 
times short enough or, as is equivalent, sufficiently large 


rates of shear to detect such effects. Under the con 
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ditions of Reiner’s experiment, the separation of the 
dises would have to be about 4 & 1077 em., which is 
about 1/25th of the mean free path, for the rate of shear 
to be comparable with the reciprocal of the relaxation 
time. 

The non-Newtonian properties can be 
according to the kinetic theory of gases (see, for ex- 
ample, Chapman and Cowling,® page 265), as an effect 
which is described by the second approximation to the 
relationship between the stress and rate of strain 
On substituting into the expression for the 


calculated, 


tensors. 
non-Newtonian properties the velocity field given by 
the low Reynolds Number solution of the Navier 
Stokes equations (1.e., Stewartson’s solution) and 
neglecting the thermal stresses (which can be shown to 
be small), it is found that the non-Newtonian properties 
of the air are equivalent to a body force 


(u2w?/2pohy”) [(@s — @2)r, O, @2(3r + 23) 

per unit volume, where cylindrical polar coordinates 
@, are numbers related to the 
to Chapinan 


(r, 0, z) are used, and @., 
molecular structure of air. According 
and Cowling, it is not unreasonable to take @ = 2, 
a = 8. It can then be shown that this body force 
gives rise to a suction at the center of amount 3y7w"a? 
2pohy?, which is about 0.02 em. of water when the gap 
is as small as | micron. Nevertheless, Reiner states 
that he disagrees with Jeffreys and that he believes 
that the very large effects he observes are due to the 
existence of non-Newtonian properties in air. 

It is important to realize that, if Reiner’s interpreta 
tion of his experiment is correct, the Navier-Stokes 
equations do not describe the mechanics of airflow. 
These equations are founded on the assumption that 
there is a linear (Newtonian) relationship between the 
stress and rate of strain tensors. The non-Newtonian 
effects postulated by Reiner to explain his experimental 
results would, if they exist, be important in other aero 
dynainical phenomena such as boundary layers or 
turbulence where similar rates of shear occur. On the 
other hand, the Navier-Stokes equations have been 
used to make predictions about these phenomena which 
have been verified in nuinerous experiments; in fact, all 
existing theories of aerodynamics which take account 
of viscosity are based on them. For this reason it 
seemed to us important to find out whether Reiner’s 
experimental results can be interpreted in the light of 
the Navier-Stokes equations without assuming that 
air has the non-Newtonian properties which he postu 
lates. 

This was the primary object of the work described 
in the present paper, but the flow situation in Reiner’s 
experiment is similar to that investigated in the theory 
of hydrodynamic lubrication, and our results seem to 
have some application in that field. 

When two parallel discs, rotor and stator, are brought 
sufficiently close together, the Reynolds Number based 
on their mutual velocity and their distance apart is small 
and Stewartson’s result, according to which the suction 
is independent of Reynolds Number, should apply. 


riICAL 


"CIENCES A4UGUST, 1957 
This result was given for incompressible fluids, but } 
seems that a similar result can be given if the fluid 
compressible. The suction at the center has a differen 
value but is constant if # varies while w and a rema; 
There is therefore no possibility of explair 
N avier 


Stokes equations if his apparatus is really functionin 


constant 
ing Reiners experimental result using the 
in the 1deal manner which mathematicians postulat 
for engineering constructions—that is, as a_perfe 
plane exactly parallel to another plane and rotating 
at uniform speed about a point in the plane. 

The questions to be considered now are what sort ; 
deviations from the ideal may be expected in apparatus 
constructed by competent engineers and what effect 
would these have on the distribution of pressure be 
tween the discs. Perhaps the most obvious deviations 
from the mathematical ideal is the rise in temperatur 
arising from the large amount of work done on th 
film of air. This can be shown, however, to be far toy 
small to produce any measurable effect on the central 
pressure. Further deviations are the inevitable lac 
of geometrical and dynamical perfection in the appa 
ratus. ‘These give rise to effects which, in turn, can b 
idealized and discussed in the light of the Navier 
Stokes equations. We shall now discuss two of these 
in soine detail: 

The first (A) is the effect of a slight error in the 
perpendicularity of the plane of either the stator or 
rotor to the axis of rotation. The second (B) is the 
effect of vibration of the rotor in the direction of its 
axis. We have, in fact, analyzed more and shall giv 
the relevant results later. 


ERROR IN PERPENDICULARITY AND 


VIBRATION 


SMALL 
EFFECT O! 


EFFECT OF 


It has been seen that the effect observed by Reiner 
was very much larger than, and of opposite sign t 
that which would arise from inertia. In the first in 
stance, therefore, inertia will be neglected, and _ th 
equations of hydrodynamic lubrication which involv 
only the effect of Newtonian viscosity will be used 
see, for example, Lamb,® section 330a). The thick 
ness /} of the thin layer between rotor and stator 1 
taken as a function of polar coordinates 7, @ and of th 
time ¢. When the rotor is rotating with angular ve 
locity w, the relevant equations for the pressure p are 


(h® 124) (Op Or) + ah 0 


(h® 12u) (Op /rd0) + ah — (1/2)hor = 0 (02 
Here wu is the viscosity, w and v are the radial and ta! 
gential components of the velocity, and the bars indi 
cate that the mean value through the thickness of th 
These equations are a 


OO 
s 


layer has been taken. 
approximation provided that (grad h)|, which is ¢ 
order ft a where a is the radius of the dises, is small 
compared with 1 and the pressure gradient due to m 
ertia, which is of order (3/10) pw*r, is small compared 


with that due to the viscous stress which is of order 


' 


0 0 


This 


pend 


00 


wher 


Both 


ill 


EPFECTS OF CORPRESSIBILI 
; 
ds but it . 5) : 
_— yor h® i.e., provided that pwh*/20u is small com 
he fluid js ee a 
— pared WiUIl I. 
1 differen ; : ae 
ere! i Che equation of continuity 1s 
a remaiy 
f ¢ < plait BO re) hii p (O O68) (hip) + r(O/Ot) (ph 0 (3 
a : nd the density p will be taken either as constant if 
: 


nctionyy : ‘ . p 
TIM the fluid 1s incompressible or as proportional to the 


OStulat " ° , _— 
| 9: pressure if it is compressible so that p/ po p/po. The 


a perfe : ‘ : : , 
ae iustification for this assumption is that the layer is so 
rotating ; oe a tia , : sae ven 
thin that conductivity prevents any appreciable vari 
ition in temperature due to adiabatic compression 

it sort 9 . ‘i000! . , ; : 
from occurring. Eliminating “ and @ between Eqs. 


\Ppparatus 29 


iat effect 


and (3), the equation for p is 


(0/086) [( ph® ur) (Op/ O86) | 
12(0/Ot 


(Op Or) } 
6w(0/O0) (prh) + 


ssure be 0 Or) |(prn’/p 
] 
eViations (prh) (4) 


perature ; . . Ke 
[his may be expressed nondimensionally by writing 


on. the 

e far to } 7/a, t tw p (p Po) pr } ie 

e central bh’ =h/ite, p po(l + p’ ( 5 

me ae kg t) then becomes, taking the viscosity as ime 

le appa pendent of the density, 

1, can be ; 

Navier- | (0 Or’) [r’h’ (1 + p’) (Op’/Or’)] + (0/08) [(h" 7) X 

of these 1 +p’) (Op’ /08@) Ar’(0/00) [h’(1 + p’)| 4 
2dr’(0/Or’) [h'(1 + p’ (6 

r in the where A Gwud*/ pyto”, and the condition to be satis- 

wapteg | fied at the edge r’ lis p’ = 0. Eq. (6) can be used 

.) 1s the , ‘ 4 ; Q te . 

aes ) represent either of the two deviations from the 


inathematical ideal. 


all give ‘ . . : 
For case (A), the error in perpendicularity, we take 


h’ 1 + er’ cos 0 (7) 
TY AND where ehy/a is the angle between the plates. For case 
B), vibration parallel to the axis, we take 
Reiner h’ | + bcos t’ (S 
—_ sai where bho is the amplitude of the vibration, ¢” nt 
first in nt’ and n/2zm is the frequency of the oscillation. 
nd the Both e and 6 are numbers which range from 0 when the 
involve plates are parallel and free from vibration to 1 when 
ve use’ | they touch. 
thick For case (A\, Eq. (6) becomes 
ator 1 
of th 0 or’) {r’(1 + er’ cos 6)* (1 + p’) (Op’/dr’); 
lar ve O OO 1 + er’ cos 0)*/r’| (1 + p’) (Op’ od ‘ 
» i \r'(2/08) {(1 + p’) (1 + er’ cosa} (9 
id, for case (B 
0/Or’) jr’(1 + p’) (Op’/or’)} [N/(1 
; CO oot"); + p) (1 b cos t”); 10) 
id tat vhere A 12nua h 
s indi 
of th 
700 INCOMPRESSIBLE FLUID 
is ( ; : , 
al lt the fluid is incompressible, p disappears from Eq. 
‘et +), and Eqs. (9) and (10) are changed only by the 
\pared mussion of the factors (1 + p’). In case (A) it will 
noticed that, since the boundary condition is p’ = 0 


order : : 
it) |, p’ must be antisymmetrical about the line 


T 


y 
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6 0. Thus the value of p’ at the point (r’, 6) be 
comes —p’ at (r’, —@). The mean value of p’ over a 


circle at any given radius is therefore zero. 


Similar considerations apply to case (B), for, taking 


of p’ at a fixed point, the right-hand 


(10) can be integrated when the factor 


the mean value 
side of Eq. 
(1 + p’) is omitted, and its mean value over a period 
is therefore zero. 

Thus, it will be seen that, if the fluid between the 
plates is incompressible, neither of the deviations (A 
or (B) from the mathematical ideal could account for 
the rise in pressure at the center which Reiner observed 
Nonzero pressures at the center will arise only if the 
lack of geometrical perfection is asymmetrical about 
every diameter—for example, if there is a spiral mark 
ing or elevation on one of the plates. But an exami 
nation of the equation for an incompressible fluid shows 
that the excess pressure then depends on the direction 
of rotation, and a pressure rise is altered to a pressure 
drop of equal magnitude when the direction of rotation 
of the motor is reversed. It follows that a mean pres 
sure rise independent of the direction of rotation and 
different from zero cannot occur if the fluid is incom 
pressible. 

COMPRESSIBLE FLUID 

It is commonly thought that the effects of com 
pressibility are important only when Mach Numbers 
are high. In Reiner’s experiment where the speed of 
rotation of the rotor was, say, 9,000 r.p.m. and the 
diameter of the discs 6.7 cm., the Mach Number of the 
fluid was less than 0.05. At first sight, therefore, one 
might not expect compressibility to alter the conclusion 
just reached for incompressible fluids that the small 
errors of construction just considered cannot give rise 
to the central pressure observed by Reiner. On the 
other hand, it has been found in the study of the me 
chanics of air-lubricated bearings that compressibility 
may have importance in cases where fluid is dragged 
into a narrow space and that the flow of a compressibl 
fluid may be very different from that of an incompres 
sible one.’ 

Unfortunately, Eqs. (9) and (10) are nonlinear and 
their complete solution would be difficult, but, if ¢ is 
assumed small in Eq. (9) or 6 in Eq. (10), analytical 
solutions of the linearized equation can be obtained 
It is also possible to linearize these equations in an 
other way. If we had considered a compressible gas 
for which up was independent of density instead ot 
would have disappeared from 
(9), and the resulting linear 


u, the factor (1 + )p’ 
the left-hand s.de of Eq. 
equation could be solved numerically using the relax 
ation technique even if e were not small. Since the 
greatest effect of compressibility seems to arise from 
the factor (1 + p’) on the right-hand side of Eq. (9 
(this is indicated by the analytical solution obtained 
below for small ¢), it is useful to have such a solution 
and Miss G. Vaisey has given the one represented in 


Fig. 1(b) for the case when the error in perpendicularity 
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Fic. 1. Pressure between dise rotating about its axis of sym 
metry and a fixed disc set at a small angle toit. Left: (a) Con 
tours for incompressible fluid with A = 4, e = 0.8, atmosphere = 
1,000. Right: (b) Contours for compressible fluid under the 
same conditions. (Relaxation solutions by G. Vaisey 


of the stator ise = 0.8, \ = 4. In Fig. 1 the contours 
represent values of 1,000 p’ 

For comparison, the contours of p’ are also shown in 
Fig. l(a) for the case when the fluid is incompressible 
(this Figure is also due to Miss Vaisey), but otherwise 
the conditions are the same as for the compressible 
case. It will be seen that the effect of compressibility 
is very large. The part of the field where there is suc 
tion which covered half the field in the incompressible 
case has contracted, and the maximum suction is now 
only about one fifth of the maximum excess pressure. 
The maximum pressure is in fact 1.6 atmospheres while 
the minimum is 0.9 atmospheres. It will be noticed 
that the effect of compressibility is to raise the average 
pressure at all radii including that at the center. The 
value § = +4 was chosen for numerical computation 
because it corresponds with that of Reiner’s rotor 
rotating at 9,000 r.p.m. when /iy) = 16.7 microns —1.e., 
near the distance at which Reiner’s effect appeared. 


ANALYTICAL SOLUTIONS 


When e or 6 are small, analytical solutions of Eqs. 
(9) and (10) can be found. Eq. (10) for case (B) will 
be considered first since the analysis is simpler. 

Case (B): We expand the dimensionless pressure 
p’ as a power series in } 


p’ = prb + pb? + psbi +... (11) 


substitute into Eq. (10), and equate the coefficients of 
b and b?, obtaining 


N(Op; ot”) = 
—Nbsint” (12) 


(1 ‘r’) (0/Or’) [r’(Op, Or’) | 


(1 r’) (0/Or’) [r’(Op. Or’) | — N(Op. Ot”) 
N(0/ot") (p; cos t”) — (1/r’) (0/Or’)r’ X 
[pi(Op, Or’) + 3(O0p,/Or’) cos t”| (13) 
The boundary conditions on /; and fp, are that they 
should be zero at the edge of the discs where r’ = | and 
finite at the center where r’ = 0. Tosolve Eq. (12), put 


‘ 


pi = f(r’) cos t” + g(r’) sin t” 
and equate coefficients of cos ?” and sin ¢”, obtain 
1 r’) (d dr’) \r’(df dr’ Ng = U0 i Th 
1 r’) (d/dr’) [r'(dg ‘dr’)| + Ni \ of 
We subtract 7 times the second equation from the fir 
and write ne 
less 
x = f(r’) - ig(? Q 
giving 
1 r’) (d'dr’) [r'(dx dr’ Nix \ 
The solution of Eq. (14) which satisfies x (1) = 
x (QO) finite is (see McLachlan* D 
x(r’) = [Jo(r’% —1iN) J ay —1N 
Writing Phe 
thos 
Jo(xV —1) = berx + ibety B 
ve f 
and Jo(xV —t)| = M(x) 
the solution of Eq. (14) is 
whe 
f(r’ [AS(V N) (berV N ber r’'V N 4 
beiV N bei r'V N) | ] 
g(r M(V N)|-2? (berV N bei r'V N — | 
le ; h 
beiV N ber r'VN = 
Litt 
We note that p; = 0, where the bar denotes a me 


at a given point with respect to time. This is to 
expected since the value of /p; is clearly unaltered if th) Hen 
fluid is taken as incompressible and also because t 
mean pressure must be independent of the sign oi 


since a change of sign of b corresponds to the additi 


of « tot”. We are, of course, interested primarily «| B 
the mean pressure since this is what is observed ¢ 
perimentally; so, taking the mean of Ea. (15), we hi 
el 
(1/r’) (d/dr’) [r’(dp./dr’)| = ! 
—(1 r’) (d/dr’)r'(d dr’) |(1 2)pi? 4+ 3p, cos 
Ww 
This equation can be integrated immediately to give — tye, 
ps = —(3/2)f(r’) — C4) [POD + g(r) radit 
where f(r’) and g(r’) are given by Eq. (15). 
Eq. (16) gives the distribution of mean pressure acr' 
the discs due to vibration. For values of .V sm Toc: 
compared with 3, it is found that a first approximatt!| grate 
to Eq. (16) is 
pb. = (N 27) (7 — Ss? i 
For .V large, p, ~ 5, 4 provided that (1 — r’)A = Subs: 
that is, the pressure between the discs is constant a! 
independent of the frequency except in the neighbor 
hood of the edges of the discs where a rapid transit 
to atmospheric pressure takes place. 
Inspection of the equation for p; shows that /; = 
so that the error is using Eq. (16) is of order 5°. radiu 
particular, the mean pressure p’ at the center 1s Ch 


EFFECTS OF COMPRESSIBILITY 


) 5 4 — lberV N/M?(v N)| —_ 
obtainigy 1 4.M2(V N)|! + O(b4) (18) 
0 ) This expression is always positive so that this effect 
on ) of vibration is to raise the pressure at the center. 
Case (.1): Here the plates are inclined at a small 


ym. the firel : : : 
“7 angle, and Eq. (9) applies. We expand the dimension- 


less pressure as a power series in e, substitute into Eq. 
9), and equate the coefficients of e and e*, obtaining 


r') (0/0r’) Ir'(Opy or’)| + 


\ | r’?) (O*p, 06") — A(Op; O8 —\rsiné@ (19 
r') (0/Or’) [r’(Op. Or’)| + (1/r"?) (07 p2/087) — 
x (= MOP: 06) = A(O O68) (r’py cos 8) — (1/r’*) X 
0 O8 [pi(Opi og) + 3r’ cos O(Op, O68) | — (j/r’ y 4 
0 Or’) [r’pi(Op, Or’) + 3r’? cos 6(Op; Or’ (?0 
[he boundary conditions on /P; and p, are identical with 
| those for case (B). 
By a method very similar to that employed above, 
we find that 
pi = lir’) cos 6 + B(r’) sin 0 
where C(7’} = 1 — 1B satisfies the equation 
l/r d ‘dr’ lr’(dC dr’ — 
C r’?) — 1AC = Kir’ (21 
\ | Phe solution of Eq. (21) which fits the boundary con- 
ditions 1s 
tes a me 
Cl | J \ 1X) Jy(V —in } 


is is to 
‘ered if tl Hence, 


ecause t 


: ; MV X (beryV X beryr’V X 4 
sign ol ] 
ie additi bei,V Xd benr’V dX) — 3 99 
rimarily B MiuV AX ber, V X betyr’V X — | 
served € bei, X he nr'~ d 
), we hi 
where J;(x¥ i beryxn + thenx and A(x 
Tey Vv ) 
OP cos 
We again note that p; = O, where the bar denotes 
- to giv . : _ 
to give | the mean value of the pressure around a circle of given 
’ radius —1.¢ 
ra 


” 


py = (1/2m) p, do 
sure acros — 


r N sm 


‘oximatt 


lo calculate p., we take the mean of Eq. (20) and inte 
grate to give 


*!] 


b —3 r’ cos W(Op, Or’ )dr’ — (1 /2)py 
. oe Substituting for pi, we then obtain 
istant al >] 
neighbor oe 'e r'(d.l dr')dr’ — (e?/4) (4? + B?) 4 
transil ; 
O(e*) (23 

—it p3 = = () } R - . , . “rel : £. 

ae since the mean pressure over a circle of fixed 
> >*. di : ° ° ° e 
er radius must be independent of the sign of e. 


- The radial distribution of the mean pressure as 
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—— A=100 \ \ 
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Ya 
Fic. 2. Radial distribution of mean pressure between rotating 

and fixed inclined discs. Curve (1) gives the values of 256p’ 


poe*d* valid for \ small compared with 16. Curves (2) and (3 
give p'/poe? for X = 16 and A = 100, respectively 


given by expression (23) is shown in Fig. 2. Curve 
(1) represents the values of 2°)" v*\* which is applicable 
when J is small compared with 16 and is derived from 
the first term in a series expansion in powers of \° 


namely, 


Curves (2) and (3) were calculated using the tabulated 
values of the Bessel functions of complex argument,’ 
and they represent the values of p’ e* for A 16 and 
4 = 100, respectively. Curve (1) gives, in fact, the 
value of p’/e* for d 16 as calculated from the ap 
proximate expression (24 
\ = 16 are about one half of these, it follows that the 


Since the actual values for 


mean pressure does not increase as rapidly as the square 
of the angular velocity of the rotor, the difference be 
coming noticeable as \ approaches 16. 

The second term on the right-hand side of Eq. | 
turns out to be small compared with the first, and, since 


9» 
=~) 


it can be seen that the second term arises from the 
factor (1 + p’) on the left-hand side of Eq. (9), this 
provides some justification for the method of lineari 
zation adopted for the solution obtained using the re 
laxation technique. In fact, the values of the mean 
pressure calculated from Miss Vaisey’s numerical solu 
tion of Eq. (9) with ¢ 0.8 and \ 
able agreement with those predicted by the expression 


tf are in reason 


(23) with the same values for e and X, the latter being 
about 15 to 25 per cent lower. For instance, the 
value of 1,000 p’ at the center is calculated from Eq. 
(24) to be 33, while from the numerical solution it 1s 
38. However, the numerical solution has a large peak 
near r’ = 1, and the agreement there with the analytical 
solution is not good. The reasons for this peak will 


be discussed in the next section. 
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Fic. 3. Calculated distributions of pressure between rotating 
cylinder and plane when separated by a small gap. Solutions of 
(d/do (p/pi) = B cos? o A( py Pp) cost & 


To conclude this section, we remark that the solu 
tions obtained above of Eqs. (9) and (10) are linearly 
superposable —that is, the effect of the deviations from 
the mathematical ideal represented by cases (A) and 
(B) may be calculated separately and added to give 
the mean pressure with error of order c*b? when there 
are both vibration and an error in perpendicularity. 

OF COMPRESSIBILITY 


COMPARISON WITH THE EFFECT 


IN AIR-LUBRICATED BEARINGS 


It will be noticed in Fig. 1 that there is a high con 
centration of pressure near the point r’ l@=f 
where the plates are closest together. In this rezion 
the air is dragged by the rotor along a circular path into 


When 


is small 


a passage which first narrows and then expands. 


the plates nearly touch—i.e., when 1 — e 
the flow in the region of near contact is like that in the 
space between a cylinder of radius R = a* hy and a 
fixed plane when the cylinder rotates with surface 
velocity dw and the minimum gap is h; = /y(1 — e). If 
x is the distance from the point where the gap is a 


minimum and /,; is small compared with KR, 


h = hy + (1/2) (x?/R) (25) 
Using Sommerfeld’s variable ¢, x = V2Rh, tan ¢, 
h = h, sec® $, and the range — © < x < +2 is cov- 


l 


as 
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if th 
AL { 


ered by —(1/2)4 < @ < (1/2)zm. Assuming tl 
temperature of the film of air is constant and eqy 
to that of the plane and cylinder, the equation for t 


distribution of pressure can be reduced to 


[( py p) cos? gq | 


ddd) (p pi B cos’ @ 


where B Oudw IR hyp, 


and p; 1s the pressure at @ Q. The required sok 
tion of Eq. (26) is one for which the pressure is equ 
to p 


mines the value of 


at +(q7/2). This condition, in fact 
A which, in nondimensional fon 


represents the rate at which air passes through the gay 


When B is very small, p p; is nearly 1 and it may bk 


taken as | in the last term of Eq. (26). It is then fou 


that 1! 1B 3 and that the solution of Eq. (26) is 
(p pi) /pr = —(1/24) (2 sin 26 + sin 46) (27 

and p; = fo since sin 2¢ and sin 4¢ both vanish at ¢ = 

+(m/2). This solution is well known in the theory 


hydrodynamic lubrication when the lubricant is as 
sumed to be incompressible, and in that case it is noi 


limited to small values of B. The excess pressur 
p — po is antisymmetrical, being positive when 4, an 
therefore x, is negative and negative when wx is positiy 
Evidently the lubricant produces no resultant pressur 
perpendicular to the plane. 

may be solved numerically for any give 


1 and B 5, ha 


Eq. (27) 
value of B. Two examples, B 
been solved by Dr. J. C. Miller. 
Fig. 3, and in the case of B 
is displayed. The calculated distributions of pressur 


5.3235 were followed throug! 


for 1 5.3220 and A 5.32. 
and it was found by extrapolation that the value of 
5.3238 gives pp) = 0.412 at @ (7/2) as well as 

o -(w 2). The pressure at the center is therefor 
0.412 


In Fig. + the distribution of pressure with x for B 


pi p 2.44 atmospheres. 
is shown, and for comparison the distribution of pres 
sure on an incompressible fluid of the same viscosit 
and under the same conditions is shown. It will | 
seen that the effect of compressibility is to increase t! 
area where there is a positive excess pressure at th 
expense of the area where it 1s negative. 
Vaisey’s calculation exhibited in Fig. 1(b) can be inter 
preted as physically analogous to the hydrodynan 
regime in an air-lubricated bearing. 
The distribution of pressure for B 1 (see Fig. 3 
much nearer to the antisymmetrical distribution ‘ 
tained with incompressible fluids, so that it is qualtt 
tively true that the main change due to compressibilit 
occurs as B to 5. Ultimately, { 
large values of B the pressure distribution becom 


increases from 1 


symmetrical. The broken curve in Fig. 4 shows t 


limiting form as B tends to infinity.* 


* An analytical solution of Eq. (26) valid for large values 
has been obtained and will form the subject of a separate not 
It is found that po/pi ~ 10.5/B? This is in g 
agreement with the value of po/p; obtained by Miller with B 
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EFFECTS OF COMPRESSIBILI* 


EFFI 1 ROUGHNESS WHEN Two SURFACES SLIDI 
Past ONE ANOTHER 

It will be noticed that the situation where the cylinder 

nd plane are close together is the same as that between 

two - vy or rough surfaces when they move past one 

nother tangentially. Thus, when two such surfaces 

ve past one another so that they are nearly in con 


tact at the high spots, the effect of compressibility is 
to make the air laver between them exert a pressure 
to hold them apart. 


tend11 


NUMERICAL COMPARISON WITH REINER’S EXPERIMENT 


Popper and Reiner* report that at 7,000 r.p.m. the 
manometer connected to the hole in the center of the 
but, the 


20 microns 


stator registered a suction, when was 


gap 
the suction dis 
If a 
so that this is just the point 


closed to 0.02 mm.--1.e., 
ippeared and excess pressure began. 3.35 em., 
this corresponds to A 2.2, 
it which the effects of compressibility might be ex 


pected to begin to show themselves if there is any error 


ef manufacture of the type considered under case (A 
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In this connection 1t is worth noticing that, if it were 
argued that Reiner’s apparatus is so exceedingly well 
made that it is unlikely that any large effect could be 
produced by a manufacturing error of the kind con 
templated here, one should consider how the experi 
The 


towards the rotor, the pressure increasing as the gap 


ment was carried out stator was moved up 


closed. If there was literally no manufacturing error 


the gap would close until / 0. This is evidently not 


a realistic approach There must be some angular 


error in the placing of the stator. Suppose it is a 
The gap can then be closed until, when / aa, there 


The 
rhus, the better the apparatus 


will be contact on the edge and ¢ l. value of 
A will then be 6u / pya? 
is constructed and the smaller therefore is a, the bigger 
will be the value of \ attained before actual contact of 
the discs. For this reason and because the mean pres 
to the 


curious conclusien that, if the pressure observed by 


sure is an increasing function of A, we come 
Reiner is due to this particular error of manufacture 
the smaller this error is, the bigger is the pressure which 


can be produced by narrowing the gap between the 
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TABLE | 
Pressure at Center of Stator When Rotor Vibrates in Direction 
of the Axis of Rotation 


(1) Ao (microns) 20) 10 5 

(2) VN 2.32 4.65 9.3 
(3) p’ b? 0.82 1.04 1. 24 
(4) 6 (with 6 = 2.5 microns 0.125 0.25 0.5 
(5) p po (atmospheres ) 0.0128 0.065 0.31 
(6) p — py (cm. of water) 13 66 320 

plates. This remark must, of course, be qualified by 


the condition that the gap should be large compared 
with the mean free path. This, however, is about 0.1 
micron at atmospheric pressure, so that the fluid may 
be treated as a continuum in the range of Reiner’s 
experiment. 

Similar remarks can be made about the effect of 
vibrations which also produce a positive pressure at the 
center of the discs. Reiner appreciated the necessity 
for reducing this error and stated that the frequency 
of the natural vibrations of his apparatus was well 
above that of the rotation of his disc so that resonance 
would not be expected. On the other hand, shafts 
running in ball bearings are always liable to suffer a 
certain amount of longitudinal vibration, and, until 
that has been measured, it is not possible to say which 
of the two inevitable deviations (A) or (B) from the 
ideal produces the larger effect at the centers of the 
discs. 
or even more, may be expected in the unfiltered air in 
Never- 


In practice, dust of dimensions of | or 2 microns, 


which Reiner’s apparatus has been operated. 
theless, it is worth while to put some figures into our 
formulas to find out whether the observed pressure is 
within the range predicted by the Navier-Stokes equa- 
tions for a well-made apparatus operated by carefully 
making the plates approach without touching. 

For a numerical example for case (B), we may imagine 
a dise 7.3 cm. in diameter (approximately of the size 
used in Reiner’s apparatus) running at 7,000 r.p.m., 
and suppose that it executes forced oscillations of small 
amplitude 6 in the direction of the axis of rotation at 
frequency 7,000 cps. In Table 1 the corresponding 
values of VN (see Eq. 10) are given in line 2 for vari 
ous gaps fy. In line 3 are given the values of p’ ‘poh? 
calculated from Eq. (18) using tables of Bessel func- 
tions. To obtain values for the mean excess pressure 
pb’ py» at the center, it is necessary to assume a value for 
6. A dise running at 7,000 r.p.m. mounted on a well 
constructed ball-bearing shaft might be expected to 
amplitude at least 


Assuming 6 = 2.5 


have longitudinal oscillation of 
0.0001 in.—i.e., 2.5 microns. 
microns, the values of b = 6 hy are given in line 4, and 
the values of the excess pressure P — fp in atmospheres 
are given in line 5. The corresponding pressures in 
centimeters of water are given in line 6. 

It appears from Table 1 that at 7,000 cps a vibration 
of amplitude of only 0.0001 in. would give rise to a pres- 
sure of 7 cm. of water when the gap was closed to 0.02 
When the gap to 5 
microns, the mean pressure produced would be 0.31 at- 


mm. = 20 microns. is closed 
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TABLE 2 

Pressure at Center of Stator When Inclined at Small Angle 1 
Rotor 

(1) ho (microns 20 10 5 

(2) VA 2 } sS 

(3) p’/e* for 7’ = 0 0.21 0.37 0. 6(approy 

1) e (with a = 4 x 10 

deg ). 125 0.25 0.5 
3) p po (em. of water 3 23 150 


mospheres, and still higher pressures would be produce 
as the gap closes still further and 6 approaches the valy 
1 which it would have when fy = 2.5 microns and th 
plates touch. 

For case (A) a similar table may be constructed 
We imagine a disc 9 em. in diameter rotating at a 
angular velocity of 7,000 r.p.m., and suppose that the 


stator is inclined at a small angle a. In Table 2 th} 


corresponding values of V \ are given in line 2 for th 
various gaps. In hne 3 are given the corresponding 
values of p’ ‘pe? at the center calculated from E 


(23). We assume, for the sake of example, that a = 
0.004 The values of e = aa ‘hy are given in line 4 an 
the values of the excess pressure p — /» in centimeter 


of water in line 5. 

It appears from Table 2 that, at an angular velocit 
of 7,000 r.p.m. and an angle between the plates | 
0.004 
the separation of the centers is 20 microns and rises 


a pressure rise of 3 cm. of water occurs whe 


to 150 cm. of water when the gap closes to 5 microns 
Still higher pressures will be produced as the gap close 
further. 

It is of interest to that the 
formulas (18) and (28) lead to finite limiting values 
Thus, in case (A) when e = 


notice approximat 
the central pressure. 
and \ ~ infinity, (pb — po) / po ~ 0.75; and, in case (B 
1 and V ~ infinity, (p bo) /po ~™ 1.25. 
must be remembered, that these 

have only been proved for small values of e and 0. 


as 6 = 


however, formula 


EXPERIMENTAL RESULTS 


We have seen that the Navier-Stokes equations wi 
predict pressures of the order of magnitude observe 
by Reiner if account is taken of the variable thicknes 
of the gap between the plates. However, it is ml 
possible to determine whether the above theory is & 
sentially correct and is the explanation of the phi 
nomena observed by Reiner until the experiment 
repeated under conditions for which the gap betwee! 
the plates may be accurately measured. 

To test our conclusions about the effect of anguls 
error, we enlisted the help of Professor Hawthorne a! 
Mr. Barker of the Engineering Laboratory, Cambridg 
We first mounted two ground and lapped steel plate 
of 4-in. diameter on a high-class boring machine ‘ 
pable of running up to 2,500 r.p.m. The lower plat 


occurred we wotl 


was insulated so that if contact 
know it and so that we could measure the gap usillt 
the plates as a condenser. The stator had five | - 


mm. holes drilled along one radius so that we cou! 


In 
perim 


mum 
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e produce 


s the value | 
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EFFECTS OF 


In 
the 


measure the radial distribution of mean pressure. 


order that the pressure we measured would be 


mean overall points on a circle, we set the stator as 
perfectly as we could, the irreducible error being about 
) microns. When the machine was set in motion and 
taking very great precautions to ensure the best pos 
sible degree of parallelism, we observed pressures up 
to 9 em. of water when the machine was rotating at 
2 OOO r.p.m but the excess pressure was greatest at 


bout half radius and nearly zero, both at the center 
We then found that one 


ind at the outermost hole. 


of the plates was very slightly conical; we lapped it 
wain, repeated the experiment, and obtained only 
0 cm. pressure. We then set the rotor at a small 
ingle to the stator so that its edge moved up and 


down relative to the stator through a distance of 2d 


0.25 X 0.001 in.—1.e., d = 0.00012 in. = 3.2 & 10-‘cm. 
fhe pressure distributions obtained with this setting 
at 2,000 r.p.m. are shown in Fig. 5. 

When running the machine the gap gradually closed 
down, owing no doubt to the slight heating of the rotat 
ing shaft, and the curves show successive distributions 
it intervals of roughly 2 min. The initial value of h 
was 15 microns. This corresponds with e = 3.215 
0.21. No doubt the value of e rose, but unfortunately 
we were unable to measure /ty) while the machine was 
running because our electrical arrangements were up- 
set by an axial vibration which undoubtedly existed. 
We not therefore know 
sponded with any except the lowest of the four curves 
After obtaining the fourth distribution, the 


do what value of e corre 


in Fig. 5. 
central pressure rapidly rose, and before we could make 
inother set of measurements the discs touched at the 
edge, the spindle bent, and the plates were destroyed. 
The value of \ corresponding to fy = 15 microns, 

2.000 r.p.m., and a diameter of 4 in. is 2.7. The 


za with e = 


pressure at the center predicted by Eq. (25 
0.21 is then about 1 cm. of water. This is consider 


tbly less than the observed pressure of 12 cm. of water. 
However, if d remains constant while /) varies, so that 


} 
au 


is inversely proportional to /) and if \ is small 


compared with 16—as 1s the case here—then, according 
to Eq. (24), the pressure at the center is inversely pro 
portional to the sixth power of fp. The discrepancy 
could therefore be accounted for if there were an error 


ol 5 microns in the measurement of fy and the actual 


gap were 10 microns instead of the measured 15. We 
do not wish to suggest that this is the sole explanation 


ol 


the discrepancy —for instance, it will not explain 
the shape of the curves of the radial distribution of 
mean pressure which are somewhat different from those 
shown in Fig. 2 especially in the later stages of the ex 
periment —but it serves to illustrate the need for really 


‘ccurate measurements of the gap between the discs 


\PPENDIX—-FURTHER ANALYTICAL SOLUTIONS 


In an attempt to explain certain features of the ex- 
perimental observations—in particular the large maxi- 


mum which occurred at the center when the experiment 
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Fic. 5. Observed distribution of pressure when angle between 


discs was 0.0033 


with inclined plates was performed—the analytical 


solution was obtained for several other deviations 


the mathematical ideal. The problem consists 
with a given choice of the function 
and (B 


where ¢ is small compared with 1, 


from 
of solving Eq. (6 
ate, 0: 8). 
h’=1+ e(r’, 0,0’ 


and the mean pressure correct to order & 


As in cases (A , it is assumed that 
(or higher 
orders if necessary iS calculated. The analysis is essen 
tially similar to that employed before when f is re- 
solved into its Fourier components with respect to @ 
and ¢’, and an analytical expression can be obtained for 
p’ in terms of these components. In general, these 
expressions are complicated and not suited for numeri 
cal evaluation without prohibitive labor, and the gen- 
eral character of the solutions is therefore best seen by 
obtaiming the asymptotic expansion valid for large 
values of \ and the power series expansion in \* valid 
for small X. 

We first considered the case of flat dises inclined at a 
i.e., h’ 1 + e(t’)r’ This is a 
and occurs when the axis of 


variable angle cos 6. 
generalization of case (A 
rotation is not perpendicular to either the rotor or the 
stator. It was found that the solution was essentially 


the same as case (A) with e replaced by a suitable mean 
value of e(t’). 


We next considered the case of conical plates in- 


=o er’ cos 6 + 


clined at a constant angle—1.e., h’ = 1 
cr’. In this case, the effect of the conicality is to 


superpose on the pressure distribution for inclined 
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flat plates an extra pressure, which is of amount 
—4,371 e?cd* po 2°- 1,225 at the center when ) is small. 
If the rotor is convex this is a suction. 

The effect of a large central hole of radius aX in the 
stator was also worked out for cases (A) and (B 
The presence of the hole changes the boundary condi- 
tions on p’; the condition p’ finite when r’ = 0 has to 
be replaced by p’ = constant when r’ = RK. It was 
found that the effect of the hole was to reduce the cen- 
tral value of the mean pressure. 

Finally, the oscillations of a flexible plate executing 
forced vibrations were considered. Three special 


cases were dealt with: (1) hk’ = 1 + 0(1 — r’) cos 
t”; here the rim of the rotor does not move, but the 
center vibrates with dimensionless amplitude }b. (2) 
h’ = 1+ b(1 — 2r’) cos t"; here the nodal line consists 
of the points 7’ = 1/2, and the center and rim of the 
rotor vibrate with amplitude } but are z out of phase. 
(3) h’ = 1+ br’ cost”; here the center of the rotor does 


not move, and the rim vibrates with amplitude 0. 

For N small, the values of the mean pressure at the 
center for the three cases are (1) ~5.2 K 10~°b?N*po, 
(2) ~—2 X 10-4b?N2bo, and (3) ~2 X 10~*b?N*pp. 
The radial distribution of mean pressure is in all cases 
comparatively flat. These figures should be com- 
pared with case (B), in particular Eq. (17), according 
to which the mean pressure at the center for case (B), 
which is the oscillation of a rigid plate, is ~5.5 X 
10-2b2N2po. 

For N large, the asymptotic radial distributions of 


ye. (2) 


mean pressure are (1) p’/po = (1/2)b°1 — r 
p’/po = b?[(5/4) — 2r’ + 2r’?], and (3) p’/po = b?[(3/4) 
+ (r’*/2)], provided that r’ lies in the range V a 
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r’<1—N~!”, Case (1) has a large central peak and 
decreases to zero rapidly as r’ > | and, in this respect 
is similar to the later experimental curves shown ;; 
Fig. 5. Case (2) has a maximum at 7’ = O and?’ = 
with a minimum in the middle at r’ = 1.2. Cay 
(3) has a minimum at the center and rises to a maximun 
nearr’ = 1. 

It is not impossible for the experimental observatior 
shown in Fig. 5 to be explained by the existence of 
vibration of type (1), but if so this would be a most cur; 
ous phenomenon because there seems to be no obvious 
way in which such a vibration could be set up. 
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Flutter of Rectangular Simply Supported 
Panels at High Supersonic Speeds 


JOHN M. HEDGEPETH* 
Langley Aeronautical Laboratory, NACA 


SUMMARY 


The problem of panel flutter of rectangular simply supported 
tes subjected to supersonic flow over one surface is treated 
eoretically. The assumption is made, and subsequently veri 
ed, that the ‘‘static’’ approximation to the aerodynamic flutter 
forces vields flutter boundaries with satisfactory accuracy for 


Mach Numbers greater than about 2. Two panel flutter analyses 
ire performed using this static approximation in conjunction with 
thin-plate theory—one employs aerodynamic strip theory, the 
other aerodynamic surface theory. The influence of Mach 
Number, dynamic pressure, panel aspect ratio, and midplane 
stress on the panel thickness required to prevent flutter is deter 


mined for extensive ranges of these parameters 


SYMBOLS 


1 = R, 2n?3?(a/b)? 
A = R, 2n*(a/b)? (note that flutter is 
usually most critical for 7 = 1) 
B = k*? + n*x*(a/b)?*R, — n'x*(a/b)4 
B = k? + n*(a/b)?R, — n*(a/b)4 
B = (a b)*R, — (a/b)4 
C(m, @ = integral defined by Eq. (A-14 
D = plate stiffness, £h*/12(1 gp?) 
E = Young’s modulus of elasticity 
Leia <t wing ams Manes = functions used in evaluation of gen 


eralized force (see Appendix A) 
= generalized force of (rs) mode on 


(mn) mode 


WV = Mach Number 

ae = midplane stresses (positive in com 
pression) 

R, = N,a?/D 

R, = N,a?/D 

R = R,/x? 

R = R,/r? 

Sim, o = integral defined by Eq. (A-10) 

= velocity of air flow 

W = complex amplitude of panel vibration 

[see Eq. (4) ] 


dimensions of panel (see Fig. 1 
= amplitude of (mm) mode |see Eq. (19 


h = panel thickness 
k = frequency parameter, wa?V yh/D 
h = b/r2 


= lateral aerodynamic loading 
= lateral aerodynamic loading due to 
(mn) mode 
ye eS = integers 
= dynamic pressure, pl’?/2 
= time 
= lateral deflection of panel 
Y, ¥, 2 = Cartesian coordinates (see Fig. 1) 


Presented at the Aer velasticity—I Session, Twenty-Fifth 
\nnual Meeting, IAS, New York, January 28-31, 1957 
* Head, Structural Mechanics Branch 


= mass density of plate 


€ = mass ratio, pa/ yh 

¢ = perturbation velocity potential 

» = dynamic pressure parameter, 2ga*/8D 

= critical value of 

7 = transtability value of \ 

uu = Poisson’s ratio 

p = mass density of air 

o = (m/n) (Bb/a) 

ww = circular frequency 

I(y = step function, equal to zero for 0 and equal to 
unity for y > 0 


INTRODUCTION 


T= AEROELASTIC INSTABILITY of aircraft skin panels 
has been the subject of a number of theoretical 
investigations.'~*! Most of these papers!—'* have dealt 
in some detail with the problem of a one-dimensional 
panel in a two-dimensional flow and have been attempts 
either to establish approximate stability boundaries for 
use in the design of flat or slightly curved panels or to 
investigate some of the apparent peculiarities of the 
panel-flutter phenomenon. The aeroelastic behavior 
of thin cylinders has also been considered,'” '° although 
in a somewhat preliminary fashion. In addition, the 
flutter of two-dimensional flat panels in a three-dimen- 
sional flow has been treated by Eisley'’ and Luke, St. 
John, and Gross.*° *! It is to the latter subject that 
the present paper is directed. 

The work both of Eisley and of Luke, St. John, and 
Gross concerns the flutter of rectangular simply sup- 
ported plates subjected to supersonic flow over one 
side; both analyses employ modal solutions and present 
numerical results for a variety of panel aspect ratios 
tor a limited range of supersonic Mach Numbers. No 
attention is paid, however, either to the higher Mach 
Numbers (greater than 2) or to the effects of midplane 
tension and compression stresses on flutter. In addi- 
tion, neither analysis is completely satisfactory for an 
isolated panel; Eisley used linearized two-dimensional 
unsteady flow theory and Luke, ef a/., while employing 
three-dimensional theory, actually solved the problem 
of an infinite-span panel separated into rectangular 
bays by equally spaced rigid streamwise stiffeners. 

In the present paper we consider again the flat rec- 
tangular simply supported plate. The important de- 
parture from the previous analyses is in the use of so- 
called ‘“‘static’’ air forces,’ for solving the flutter prob- 
lem. Experience has shown that this static approxi- 
mation yields flutter boundaries that exhibit excellent 
agreement with those yielded by more refined (with 


—— 
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Rectangular panel and coordinate system 


Fic. 1. 


regard to unsteady effects) aerodynamic theories, 
provided that the Mach Number is sufficiently high.* 
The use of the static air forces has the obvious virtue 
of simplicity and allows the inclusion of three-dimen- 
sional aerodynamic forces in a relatively simple manner. 
In addition, large ranges of the pertinent parameters 
of the problem can be readily investigated. Also, stud- 
ies of the effects on the flutter boundaries of the mid- 
plane stresses that might arise as a result of loading 


and aerodynamic heating are feasible. 


STATEMENT OF PROBLEM 


The configuration to be analyzed herein is shown in 
Fig. 1. It consists of a simply supported flat plate of 
uniform thickness / mounted in a rigid wall with air 
flowing over the top at a Mach Number J/. The plate 
has a length a and a width 6 and is subjected to constant 
midplane force intensities V, and N, (positive in com- 
pression). 

In the following analysis, small-deflection thin-plate 
theory is assumed to apply. The appropriate equilib- 
rium equation and boundary conditions are 


DV*w + N,wWr: + NiWy + yhwr I(x, y, t) 
w(x, 0 t) = w(x, 6, t) = w0, y, t) 
w(a,y,t) =O)" (1) 


Welz; 0, 2) Wyy(x, B, t) w,,(0, ¥, ¢) 


w,(a, y, t) = 0 


where w(x, y, ¢) is the lateral deflection of the plate 
and /(x, y, ¢) is the lateral load per unit area due to 
aerodynamic pressures; D is the plate stiffness 
E:h® 12(1 — w*) and y is the density of the plate material. 

As was mentioned in the Introduction, we will assume 
that the air forces yielded by linearized static aerody- 
namic theory give an adequate approximation for this 
problem. In this theory, the loading at each instant 
is taken to be the loading which would result from flow 
over a stationary surface with a shape the same as the 
shape of the deflected plate at that instant. Thus, no 
time-dependent effects are considered in either the 
potential equation, the boundary conditions, or the 


* Comparisons presented in a subsequent section show that 
the lower limit on Mach Number is somewhere between +/ 2 and 
9 
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pressure-potential relations. Note that this 
(on the side of simplicity) from the so-called “quasi 
steady” theory, in which time derivatives are includ 
in the boundary conditions and _ pressure-potent; 
relations, and from the first-order theory,?? in why 
time-dependent effects are included approximately ; 
all three governing relations. 

In the following sections a closed-form solution of E 
(1) 1s first carried out for two-dimensional aerodynan 
theory (strip theory). The same problem is _ th 
solved by the Galerkin procedure in order to invest 
Finalh 
a modal solution of Eq. (1) for three-dimensional theor 


gate the convergence of this modal approach. 
(surface theory) is presented. 


ANALYSIS 


Aerodynamic Strip Theory 
For strip theory, the lateral loading is given by th 


simple Ackeret value 


Hx; 951 — (2q/B)w,(x, y, t 


where g is the dynamic pressure pl’*/2, and 3 


VM? — 1. Eq. (1) then becomes 


DVtw + NiwWr, + NywWy + (2¢/B)w, + yhwu = 0 


@ 
For the simply supported plate the solution of Eq. (2 
can be written in the form 
wW Re[W,,(x) sin (ny /b)e™ 
The frequency w is, in general, complex; however, 
direct our attention primarily to real values of o fa 
which the motion is harmonic. 


Substituting Eq. (4) into Eq. (3) and the associat 
boundary conditions and nondimensionalizing vields 


W, ° — 2n?r2(a/b)?W,,” + nirit(a/b)tw, + 


R,W,,"” — Ryn?x?(a/b)?W, + AW,’ — Rk? W, = 


W,,(0) W,,”(0) W,,(1) W,"(1 0) 


where the primes denote differentiation with respé 


tox, = x/a, and 
k= Now /D, Kt, N,a*/D, | 
r 2qa*/BD,  k? (yha*'D) w*' 


Eq. (5) and the boundary conditions constitute 4! 
eigenvalue problem where k? can be considered to! 
the eigenvalue. We seek, in fact, the values of A, f 
and R, for which one of the eigenvalues k* becomes con 
plex and, hence, the panel becomes unstable [see E 
(4) }. 


Eq. (5) can be rewritten 


w,'’ + AW,” + rW,’ — BW, = 0 
where 


A R, — 2n°x?(a/b)? | R 
B = k? + n°nr2(a/b)?R, — n'ri(a/b)*$ 


Note that B now takes on the role of the eigenvalue 
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ic. 2. Influence of airflow on panel frequencies 
Exact Solution—The solution of Eq. (7) can be ob- 


tained directly from the analysis'* of the limiting case 
of the infinite aspect ratio simply supported plate. In 
the present notation, the condition for a nontrivial 


solution of Eq. (7) becomes 


e + 6 + 4a°(6? — e€*)] sin 6 sinh e 
Sa’«d(cosh € cos 6 — cosh 2a (Y 
where 
6 V (A (4a , 2 (A 2) 
€ V (\ 4a) — a® — (A/2) (10) 
B (A° 16a" = tla? + (A 4) }* 


The eigenvalue B can be computed for fixed values 
of \ and A by varying a until Eq. (9) is satisfied and 
then using the last of Eqs. (10) to determine B. The 
resulting solutions can be shown in the form of plots 


of \ versus B for various values of A. An example 


plot for A 0 is indicated by the solid line in Fig. 2. 
The abscissa is \B/z instead of B in order to allow 
the inclusion of large ranges of B conveniently.) 

It can be seen that for a fixed value of \ there exists 
an infinite discrete set of values of B. For 0, the 
values of B can be shown to correspond to the natural 
Increasing the 
the 


frequencies of the plate in a vacuum. 
dynamic-pressure parameter \ merely changes 
values of B (and, hence, the natural frequencies) 
smoothly until a critical value (denoted A,, in Fig. 2) 
At this value, two of the roots for B be- 
For \ < X.,, the values of the frequencies 
(provided the panel is unbuckled for 
\ = 0) and, thus, the panel is stable. For \ > X,,, two 
of the values of B become complex; hence, in view of 
Eqs. (4), (6), and (8), the panel has at least one un 
Therefore, the value X, 
It should be noted that 


is reached. 
come equal. 


reiain real 


stable mode of oscillation. 


defines the flutter boundary. 
this value is characterized by the relation 


O\ OB OA ‘OR? = 0 (11 


Where the partial derivatives are performed while 


keeping A constant. 
rhe variation of \,, with A has been calculated and is 


shown by the solid line in Fig. 3. The abscissa here is 


— 2n*(a/b)? (12) 


| A/r? R, 


bIM PL 
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Xt. R,/x? 


where 


Note that, for 4 < 5, \., decreases monotonically with 
A; thus, for fixed values of R, and a 6, the lowest value 
OF Az, The resulting 


shape of the flutter mode in the y-direction would be a 


would be achieved for a 


half sine wave. 

Although a full discussion of the flutter boundary in 
Fig. 3 is included in a later section of the paper, a men 
tion of several salient points is warranted here. First, 
the flutter speed is independent of R,. This result, 
which comes about from the lack of coupling between 
the various modes in the y-direction, is of interest in 
connection with wing panels for which the greatest 
stress is usually in the spanwise (y) direction. Second, 
the range of validity of the results in Fig. 3 is limited 
by the buckling characteristics of the plate. For in 
stance, a square plate with no stress in the y-direction 
(R, 0) would buckle (for A 0) if R, t; for this 
reason, the curve in Fig. 3 would be valid only for 
A <2.* Other values of a/b and R, would yield other 
ranges of validity. 
parameter is finite when the panel is on the verge of 
buckling, even though the of the 
panel (as measured by the lowest in-vacuo natural 


Finally, the critical dy 


Third, the critical dynamic-pressure 
“effective stiffness’ 


frequency) approaches zero. 
namic-pressure parameter goes to zero at A equal to 


5. The need of large panel thicknesses if A ~ 5 is 
obvious. 

Galerkin Solution--In order to investigate the con 
vergence of the modal procedure, a Galerkin solution 
of Eq. (7) will be performed. 

Let 


Qo Cmn SIN MIX, 


Substituting into Eq. (7), multiplying by sin rmx, and 
integrating yields the following set of equations for the 


coefficients C,,,- 
(m4 — m?A 


(m a ee 13) 


* Since the buckling load is affected by the air forces, this in- 
equality is not strictly correct. This phenomenon ts discussed 


in Appendix B 
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where term solution yields good results over the entire rang, 
5 ' eens . A of A considered.* 
B B/ x‘ (k2/m') + n*(a,b)?R, — n*(a/b)3 


Riakey-2 (4/mr) [rm/(r? — m?*) |; 

n s, m-+rodd)_ (14) Aerodynamic Surface Theory 
Linn, + 0; n= 5s, m+ reven The results in the preceding section were obtain Z 
Linn, rs = OF n# on the basis of aerodynamic strip theory. The jy) 5°" 
In the first of the latter expressions, R R./x2. The fluence of three-dimensional aerodynamic effects on th. ©" 

‘ « P, Ns Sy 7) 7) . << P ? sé : ; 9 b 
seemingly overcomplicated form for the generalized flutter boundary is of interest. In this section, t} 
aerodvnamic-force coefficient F,,,.,, is introduced at flutter of the rectangular panel in a three-dimensiong) “ 
this point to agree with the notation in the next section. flow is analyzed by means of the Galerkin solution, es 

— e e * Ad ay -p de > ‘ > ‘ $ > 7 CWO 
For a nontrivial solution, the determinant of the Air Forces ~The lateral loading due to the supe | 
. c . . DOSS 
coefficients of ¢»», must equal zero. If 7 2, for in- sonic flow above the panel can be found in a straight) 
stance : forward manner by means of source-superpositic; Valu 
. . ? . ce hel [a . - py 
: is techniques. Thus, if the deflection of the panel Ps 
= io — (8/34 94 valu 
A B (S/o) _ () w(x, y, t), the lateral load /(x, y, ¢) is*4 - 
(8/3r')X 16—4A —B taine 
I(x, y, t) pU(d¢g/Ox) (x, y, b) lj It 
Upon expanding and solving for \ we get aaa 
where sets 
d = (341/8) V(A + B — 1) (16 — 44 — B) (15) : . 9s 
| “ U ¢f (Ow O&) (E, n, t) dE dn __— mids) 
Sy ae a sill as “ ” . oe ee - -- 7 ») QAd, V, — 1\ 
1 his relation is illustrated by the dashed line in Fig. 2 JJ V(x — 8)? — By — 9)? follov 
for A = 0. The two-term approximation to the critical : 
value of \ can be obtained by finding the maximum; The integration is carried out over the forward Ma 
the result is cone from the pe int (x, y). 
. ° . . | 
eo . z In the modal solution the deflection is expressed | 
Ner = (9x'/16)|5 — A (16) 
the series 
This approximation and the three-term and four-term . : ! 
PI ‘ i ; ; ene) ig w R, Rude Crs sin (rmx /a) sin (sry/b)e™ 1) Fy [me 
approximations (j7 = 3 and 4) are compared in Fig. 3 rs 
with the exact result. It can be seen that the two- ; by ak 
é i : ; in the region0 <x <a,0 <y <b. The correspon 
term approximation consistently underestimates X,,, ‘ . ; ; F 
; : + : ‘ ing loading is O 
the error increasing as A becomes large negatively. § cnet 
. Ac; CO 
The three-terim approximation is satisfactory for small l(ix, y, t) = R, [ De decradea(x ye 204 deter 
values of A but is again unsatisfactory for the larger whicl 
negative values of A. On the other hand, the four- where St 
R,, R 
ro (*® Pyt(/8(e-9 cos (rrt/a) sin (Ern/b) becor 
ls(x, y) = —2q 7 [1(n) — 1(m — b)| dn dé (2)\) real. 
aox Jo Jy-c/a)(x-p V(x — &)? — By — n)? pa 
Note that the step-function expression lor ) 
In 
I(n) — 1(n — 6) = 0; 1 < O 
matic 
L; O< 97 <6 
U] 
each 
0; n>bd 
expan 
excludes the regions outside of the plate from the integration. sultin 
Modal Solution—Substituting Eqs. (19) and (20) into the equilibrium equation, Eq. (1), gives Phe « 
proces 
— rx? s2e?\2 pkg? gtx? : rmx. osmy j i 
pe BE | > ( Sk =e = i vhs | Crs Sin sin —~ — Cysl;5(x, y) 0 
os a” b a b | a b f virtuc 
, ee = denth 
If this equation is multiplied by sin (mmx a) cin (my 6) and integrated, the recult is dine 
{[m? + n%(a/b)?]? — mR, — n*(a/b)*Ry — BY can = (/4*) DD Lenn, slr » | couph 
) ol me 
where k = k/r’ and result 
¢ ' be Cor 
2 Br (° (7 |. max . ney jO (* (2*°/"%-9 cos (rré/a) sin (snd) ). TI 
smn, 7s = sin sin ) , [1(n) — 1(n — b) |dn dé pax l¢ 
rab J0 0 a b Ox JO J y—(1/B) (x—&) Vv (x — £)? —_ B>(y a n)* ) sibilit 
: ‘ ; - : paran 
is the generalized force of the (rs) mode on the (mn) mode. Integration by parts yields 
* It should be mentioned that the four-term solution would also be inaccurate for very large negative values of A. In this casé,' 
as it 


behavior of the plate approaches that of a membrane for which the Galerkin solution is known to be inapplicable.** 





ntire rang, 


re obtain 
The j 
fects on tt 
ection, th; 
limension 
ution. 
the super 
a straight 
perpositior 


e panel | 


vard Mact 


‘pressed 


twl 1 
e ] 


orrespon 


l, 
n dé pax 


lis case, ! 


| 


| 


FLUTTER OF 


RECTANGULAR 


Ss 


{3mr [ in [ [ 8) (*~®) cos (marx a) sin (nwy/b) cos (rrt/a) sin (san b) 
ab JOZO0I0e 1/8 x 3 i 


V/ 


Before continuing, it is necessary to determine these 
generalized forces. As is usual for supersonic aero- 
dynamic problems, the generalized forces can be shown 
to be functions of the modified aspect ratio parameter 
For Bb a > 1 

in terms of Bessel and Struve functions and 


3b ‘a. , the reduction of the expression 
for Lins 
two readily evaluated single numerical integrations is 
possible; this reduction is carried out in Appendix A. 
Values of L,,, i. & 

2 3, 4 are given in Table 1. 


values for the two-dimensional case (8b a x 


for Bb a f and m, #, 7, s 
Also given are the 
ob- 
tained from Eq. (15). 

It should be noted that L,,, 


O when » + s 1s 


odd. Thus, flutter modes are either symmetric (7 
|,3,...) or antisymmetric (7 2, 4,...) about the 
midspan. Also, the generalized force possesses the 
following symmetry properties: 

| bs 

os - m + reven 
‘im? + n2(a/b)2]? — m2? R, — n2(a/b)? Ry, — R*} Cnn = 


of the 
When expanded, this 


For a nontrivial solution, the determinant 


coefficients of c,,, must be zero. 


| determinant yields an eighth-order equation in k?® in 


as it was for the strip-theory solution. 


which the coefficients are dependent on \, R,, R,, and 
bh. As before, we seek the lowest value of \ for fixed 
R,, R 


become complex. 


and ab for which two of the eight roots for k? 
(For \ = 
The calculations have been carried out for a/b 

Values 


0, all roots are obviously 
real. 
| and for a number of values of R, and R,. 
Table 2. 
1, 2, and 


for \,, were obtained as shown in 

In the calculations for 8b a $!, an auto- 
matic high-speed digital computer was used. For 
each set of values for R, and R, the determinant was 
The 


sulting eighth-order equation in k? was then solved. 


expanded for various selected values of X. re- 
The critical value of \ was found by an interpolation 
procedure. 

A significant outcome of these calculations is the 
on the value of R,. Evi- 
dently the coupling between the various modes in the 
= jf. the 
coupling is so slight, it is improbable that the addition 


virtual independence of d.¢; 


y-direction is slight, even for Bb/a Since 
ol more modes in the y-direction would change the 
results appreciably. Thus, the tabulated values can 
be considered to be essentially correct. 

The lack of dependence on R, also suggests the pos- 
sibility that the value of A., is dependent only on the 
parameter 


A= — 2(a/b)* 


R, 


If this were 


(xX 


— £2) 
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— py = 9)" 
[1(n) — 1(n — b)| dn dit dx dy (23) 
a ra) ee m + r odd 


For Bb a 


duction can be made. 


1, a similar, but more complicated, re 
No results have been obtained 
for this case, however, primarily because they would 
be useful only for very low aspect ratio panels. For 
example, values of 8b a < 1 would correspond to Mach 
Numbers less than vy 2 for a panel of unit aspect ratio 
But, as will be demonstrated, the static aerodynamic 
theory cannot be expected to yield good results for this 
region. 

Calculations for an Example Panel—-In view of the 
results obtained for the modal solution for strip theory, 
it is evident that the use of four modes in the x-direc 
tion should give satisfactory accuracy. In the y-direc 
tion, it is convenient to include two modes-— b 
the 


in order to include 


since this constituted the most critical value for 
strip-theory solution; and » 3, 
effects of spanwise coupling. 

Eq. (22) then becomes 


—~ —_ -» m |, 2, 3, 4 
nr) > 2, Lame ta *'@ 
r=1,2,3,4 1,3 n= 1,5 


the cace, the results for a b | could be extended to 
apply to other aspect ratios through the use of this 
As a check, the values of \., were calcu 
lated for a/b = 2, R, 1, R, 1/2, 
R, —3.5, R, 0. For gb/a 1, the resulting 
values of \., are 647.7 and 626.6, 1espectively, in rea 


parameter. 
0 and a b 


sonable agreement with the tabulated value for a/b 
1, R, = —2, R, = 0. 
ently be made. The variation of A., 


Thus, the extension can appar 


with 4 is shown 


in Fig. 4 for Bb/a 1, 2, 4, and © (strip theory 
the values are also given in Table 3. Note that for 
A = 5, the critical dynamic-pressure parameter is zero 


for all values of 8) a. 

The significant results of this paper are embodied 
in Table 3. From this Table, the critical values of A 
for isolated rectangular simply supported panels can 


be found in terms of A for a large range of panel aspect 
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ratio, midplane stress, and Mach Number. The 
necessary panel thickness to prevent panel flutter can 


then be calculated from the value of \.,.._ In this regard, 
it is interesting to note that the aerodynamic effect of 
finite aspect ratio on the panel thickness is very small. 
For the range considered the variation of \,, between 
Bb a x 
the corresponding variation of thickness would be less 


and 8b/a | is approximately 6 per cent; 


than 2 per cent. 
DISCUSSION 


Accuracy of Static Aerodynamic Approximation 


The results obtained by Eisley'® by using linearized 
two-dimensional unsteady flow theory provide an ex 
cellent basis for evaluating the accuracy of the static 
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aerodynamic approximation employed in this paper 


The results calculated by a two-mode Galerkin solutjo; 
for unstressed plates of aspect ratio 4, 2, 1, and 
subjected to J = 2 flow are shown in Fig. 5. 7) 


abscissa in this Figure 1s the well-known mass ratio 


which does not appear in the present analysis becay 


of the simplicity of the aerodynamic theory. Inde 
the values of \,, obtained from the static approximatj 
are the same as those in Fig. 5 for e 0. 
apparent that the static approximation is adequate 
eis small enough. 

In order to define the practical range of €, it is cor 
venient to make use of the followinz combined | 


rameter: 


TABLE 1 


Values of Generalized Force (Limn,rs) 
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r/é* 12(1 — pw?)(pc?/E)(y/ p)*\ M7/B) 

where c is the speed of sound in the air. Note that, 
for fixed Mach Number, this parameter is dependent 
only on the properties of the air and the plate material. 
For a particular altitude-material combination, there- 
fore, § would vary as e*. The variation for magne- 
sium at sea level is indicated by the dashed line in Fig. 
5. For this altitude-material combination the error 
of the static approximation can be seen to be negligible. 
For higher altitudes or denser materials the corre- 
sponding cubic line would lie above the one shown and 
the error would be less. It is reasonable to conclude, 
therefore, that the static approximation yields satis- 
factory accuracy for AJ = 2. 

Investigations of the frequency expansions of the 
flutter forces (see reference 6, for example) show that 
the accuracy of the static approximation would be even 
better for Mach Numbers greater than 2. For J/ < 2, 
however, the accuracy decreases. At J/ 4/2, for 
instance, the unsteady flow theory yields two flutter 
boundaries,® '* whereas the static theory predicts only 
one. It is probable, therefore, that a lower limit, 
somewhere between Vf2 and 2, exists on the Mach 
Number range for which the results obtained in this 


paper are applicable. 


Panel Thickness 


Two examples of the variation of the critical thick- 
ness ratio with Mach number as obtained by the pres- 
ent theory are indicated by the solid lines in Fig. 6. 
The results apply to unstressed aluminum panels at 
an altitude of 50,000 ft. with aspect ratios of one and 
infinity. As could be expected, the effect of reducing 
the aspect ratio is to reduce the required / a. 

Also shown in Fig. 6 are the results calculated from 
unsteady flow theory for the lower supersonic Mach 
Numbers by Luke and St. John.*! results, 
which are indicated by the dashed lines, merge well 
The 
meeting is nearly exact for the aspect ratio of infinity 
if account is taken of the fact that the dashed lines were 
obtained from a two-mode solution,*! whereas the solid 


These 


with the predictions of the present analysis. 


lines were obtained from a four-mode solution; the effect 
of going from four to two modes is to increase the thick 
ness by 7.5 per cent. For the case of unit aspect ratio, 
the results are again not directly comparable; in addi- 
tion to being based on fewer modes, the dashed line 
applies to a spanwise array of panels rather than to an 
The effect of this latter disparity is 
1.5, however, since strip 


isolated panel. 
probably small for JJ ~ 
theory yields the same results for both configurations 
and the aerodynamic influence of finite aspect ratio 
is rather unimportant. The correction for number 
of modes mentioned above would again produce ex- 
cellent agreement. 

Past analyses of two-dimensional panel flutter have 
indicated that the most critical condition for design 
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occurs at a value of Mach Number of about 1.5. Tp; 
is illustrated strikingly by the infinite aspect ratio cury 
in Fig. 6. The required thickness decreases rapid) 
between .\/ 
slowly. Even at J/ 
about one-fourth of that at 1/ 1.3. 
situation is apparently different for finite aspect ratio; 


1.3 and 1.5 and then increases yer 
6, the thickness ratio is on} 
However, t} 


a virtually monotonic increase of required thicknes 
with Mach Number is experienced for an aspect ratio, 


1 


one. Thus, the high Mach Number theory of th 


paper should provide useful design information fo; 


panels of moderately low aspect ratio. 

The numerical values of the critical thickness rat; 
in Fig. 6 are fairly small and would be of consequen 
only for lightly loaded portions of the structure, suc 
as fairings and control surfaces. It should be en 
phasized, however, that large thickness ratios would b 
required if the parameter A were close to 5. This sity 
ation would not be difficult to attain. For instance, 
a panel of aspect ratio equal to 1/2 were subjected to 
compressive midplane stress in the chordwise directio: 


without buckling 


only, A = 5 could be achieved 


Even for a panel of unit aspect ratio, A 5 could b 
obtained without buckling by compressing in the 


direction and pulling in the y-direction. 


that such configurations should be assiduously avoided 


in design. 


CONCLUDING REMARKS 


The panel flutter analysis presented herein has bee: 
restricted to the problem of an isolated simply su; 
ported plate of uniform thickness. The same type 
analysis can be applied, however, to other panel cot 
figurations. Clamped panels, integrally stiffened pai 
els, arrays of panels, and others can be treated by th 
modal approach based on the static aerodynamic aj 
proximation; the results should be accurate for th 
higher Mach Numbers. 


The emphasis in this paper has been on the flutter 


of flat unbuckled panels. There is a_ considerabk 
amount of interest in the flutter of panels buckled b) 
loading or heating. For this reason, the results of 
preliminary investigation of the flutter behavior 0! 
buckled included 


Appendix B. 


simply supported panels are 


It is generally agreed that flutter analyses should b 
supported by experiment. In this regard, the sma 
amount of data’ ** which is available indicates tha! 
are somewhal 


panel-flutter calculations 


conservative. Most of 
either to buckled panels or to the lower supersot! 


More experimental data, especial) 


theoretical 


these data, however, relat 
Mach Numbers. 
on unbuckled panels, are needed in order to assess U 


validity of panel-flutter analyses. 


APPENDIX A—EVALUATION OF L,,, 


92 
(23) 


The generalized force L,,,, ;; is given by Eq. 
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Lomr ~ 


iT hk 


™ 


It ¢ and 7 are replaced by 


the result 1 


: mr sas si 
of ad b /70 #0 


where the order of integration has been changed. 
after some manipulation 


When » + 


s is odd, 


( 1/8) (x-—é COS (MTX QA) 
1/8) (a 
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sin (may 6) cos (rata) sin (san /b 
¥ (* — £)° — FY — 9 
[Z(n) — 1(n b)| dn dé dx dy A-] 
Le =<) — Py n 
x=) + fy 9 
v—Uu 
y= x 
28 / 
v— u ‘ marx rn u+i 
— dy cos cos t- ax 4-2 
28 J (u 2 a a 2 


Carrying out the x- and y-integrations yields the following result 


ie 0 4-3 
When » + s 1s even, 
* ! mInrs - mz Pr 
be = a 7 (—1)" (ag — Tin — Ts + de “ 
r (ry? — m*) (s? — Hn n#s 
? ns P ; ’ m r 
(= Da — Kis A Big) 
T Ss n- Ns S§ 
A-4 
2 mr _ mez r 
: (—1)" an — Len — Jen + J 
rr? — mm n S 
l P ; m ? 
—(—1)" (— Linn — Jinn + Kirn + Linn 
~ n S 
In these expressions, 
I"m (( sin mr|(u + v)/2a| sin nr[(v — u)/2Bb| 
du dv 
aon J, V uv 
R 
m a Fs . ] ° ‘ ‘ 
1) }1 — [(v — u)/26b)} sin ma[(u + v) 2a] cos nx[(v — u) /28b 
mr du dv 
a . V uv 
R 
m ) a . ‘ ° ‘ 
: mr }1 — [(u + v)/2a]}} cos mr[(u + v) 2a] sin nmx[(v — u) 28d] : 
K, du dv A-5 
a ee V uv 
R 
— 1)" rr 
L (ma)* 
a y a 
R 
j ‘ ‘ ‘ [ ‘ 
(1 — [(a + v) 2a]j 41 — [(v — u)/26b]} cos mm [(u + v) /2a] cos nx[(v — u) 286] 
7 du dt 
V uv 
Che region of integration R is bounded by the lines Then /,,,,, for example, becomes 
u=0,u =v, and wu + v = 2a for Bb/a > 1 and by the 
additional line v — x 286 for Bb/a < 1. *r/2 Pl 
. . 7 ° . ») ZT . ~ 
Hereinafter we confine our attention to values of Im (2m n) (—1 sin may X 
° J J 
e0 a greater than or equal to one. For this range 
the expressions for J,,,, Jon Kun, and L,, can be re- sin [(mm o)f sin 6|d0 dg (A-7) 
} . . ° @ 
duced to simpler forms, as follows: 
Let where o (mn) (Bb a) (A-8) 
u = at(l1 — sin 6) \-6 
. > (A-0) = . aid , ‘ 
Uu a¢(1 + sin 6) Carrying out the ¢-integration yields 
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] a0 | ('"? sin [(mm oa) sin 6] dé or a —(a nm) [S(m, «) — S(m, —o \ 
nmr Io o + sin @ 
where 
**/2 sin [(m ao) sin 6] dA) : *7/2 sin [(mmr a |) sin 6] dé 
: . S(m, o : A-] 
Jo o — sin 6 f Jo sné-+ a 
By a similar procedure we find 
ve [o/(o? — 1)| (1/mm) [S(m, o) — S(m, —2o) | a[C(m, «) — C(m, —o 
[ao (oc? — 1 mn) [2 rll,(mr a A-] 
| —[o*/(o* — 1)] (1mm) [S(m, o) S(m, —a)| — o(mn) [C(m, «) + C(m, —a)]| - 
[o/(o7 — 1)|] (mn) [2 — rl (mr o A-] 
= [30° (07 — 1)*] (1mm) [S(m, 0) — S(m, —o)|] + oma[S(m. oc) + S(m, —o 
{o*/(o? — 1)] [C(m, o) C(m, —a)] + [o/(o? — 1)] (m/n) [C(m, o) + C(m, —o 
lo ‘(07 — 1) |mar*J, (mmr co) + [o(2 + 0°) ‘(02 — 1)?] (m/n) [2 — rH (mr ‘o)] 4 
[2a (co? — 1 mn =% A-l: 
’ **/2 cos [(mr/|o|) sin 6] — cos mr 
where C(m, o : de A-] 
J0 sin 6 + @ ‘ 


and J, and 77, are Bessel and Struve functions. 


The functions S(m, «) and C(m, a) can be evaluated accurately by means of numerical integration; the presence 


the term cos mm in Eq. (A-14) ensures well-behaved integrands. 
Table 1, Gaussian integration with ten points on the interval (0, 7 


integrals accurate to at least five significant figures. 


Except in the special case o 1 the values of L,,,,, -, can be calculated from Eqs. (A-3), (A-4), and (A-9) throug 
(A-14). Foro 1, it is convenient to go back to Eq. (A-5) and obtain expressions for Jpn, Sinny Amn, ANd Lyy di 
rectly. Thus, since m2 /a = nw 8b for this case, we can write, for example, 

sin ma[(u + v)/2a] sin nx [(v — u)/ 28d] (1 2) [cos (mru/a cos (mv ‘a 
and the necessary integrations can be performed without changing variables. The results for ¢ 1 are 
) om (m/n) [2 — wlT\(mz) | 
Tics (2/3) (m/n) + (1/3) [(mr)? n|Ho(mr) — mr?S\(mr) — (2/3) (ma n)N(mr | 
| ee 2(m/n) (—1)" — (2/3) (m/n) + (273) [(mm)?/n] Ho(mr) — (1/3) (max n) (mr A-lj 
| see (2/3) (m/n) (—-1)" + (2/15) (min) + (2/3) m2? x? Jo(mar (1/3)mnr*J\(mr) — 
(4/15) (m3x?/n) + (1/5) [(mm) /n]? Ho(mmr) + (mm’n) [(2 15)mPx? — (2.5) Aimar 


FLUTTER OF BUCKLED SIMPLY 
SUPPORTED PANELS 


APPENDIX B 


It is generally agreed that buckling has an adverse 
effect on the flutter speed of panels. Past theoretical 
work on this subject, however, has dealt with only the 
This work 


was pioneered by Isaacs,' who introduced the concept 


two-dimensional panel-flutter problem. 


of the ‘‘transtability’’ speed, a speed calculated from 
purely static considerations that constitutes an upper 
limit on the flutter speed of the buckled panel. Hayes? 
and Miles* then demonstrated that the transtability 
speed is indeed equal to the flutter speed for small per 
turbations about the static equilibrium position. 

The similarity between the governing equations of 
and three-dimensional flutter for strip 


two panel 


theory with static air forces [see Eq. (7) and refer- 
ence 12] indicates that the flutter speed of a buckled 
finite aspect ratio panel might also be obtained through 
a transtability type of approach. A preliminary in- 
vestigation, limited by the use of small-deflection plate 


theory to small buckle depths, of the static buckling 


For the calculations of the values of 1 
») 


a) 


behavior of panels is outlined here. Since some of th 


results obtained are of uncertain validity, being w 


supported by a dynamic analysis, the investigation 1 


limited to a two-mode solution for simplicity. 
The effect of the air forces on the buckling load can! 


determined by setting k equal to zero in Eq. 
Slightly rewritten, the equation becomes 


(SA/3r*)? = (A G 
where, for 7 5 


B, Ra b)? — (a b 


The buckling loads can be calculated from this equ 
An example plot of the variation of R, with 
1, &, 
would result for other values of a b and R,. 


tion. 
for a/b 0 is shown in Fig. 7. 
The transtability analysis for the two-dimension 
problem is based on a plot like that in Fig. 7. It* 
argued that a buckled panel subjected to air flow w yu 
exhibit the midplane load given by the left-hand bran¢ 


of the curve. If A is increased, the load would increas 


was used; this procedure yields values of the 
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FLUTTER OF RECTANGULAR 


maxunum of the curve is reached. Further 


until the 
increase of A would necessarily produce motion of the 
panel since there would be no stable static equilibrium 
position of the buckled panel. (Presumably, the panel 
be unstable at the higher buckling loads not con 
The 


the curve is the transtability 


would 


sidered in the two-mode solution. maximum 


value of \ given by 

luc 

The same line of reasoning can be used for the three- 
dimensional problem, provided that the buckle depth 
is kept small. For this case, the results are dependent 
on the boundary conditions on in-plane deflections 
ind stresses at the edges. Consider, for instance, a 
buckled panel with streamwise edges having a pre 
scribed stress and spanwise edges having a prescribed 
displacement. For such boundary conditions the value 
of By would be fixed and the value of A would be that 
given by Eq. (B-1). The resulting transtability value 


of \ would be found from the condition 
O\ 0A = O 


which gives 


7 (994/32 $+ 5 (B-2 


On the other hand, a panel with the edge boundary 
conditions reversed would have a fixed value of A and 


the appropriate transtability condition would be 
Or\ OB, = 0 
Phe resulting value of Av is 


Ar = (9r4 16) 5 — A (B-3 
which is the same as that obtained for the unbuckled 
panel |see Eq. (16) ]; this result can, of course, be quite 
different from that given by Eq. (B-2). 

In a practical situation, the edge boundary conditions 
would probably be different from either of the foregoing 
cases, and different values of A; would result; each 
situation requires separate treatment. 

A question might arise as to the usefulness of the 
small-deflection theory in obtaining results pertaining 
to finite aspect ratio buckled panels. It is felt, how 
ever, that such results represent correctly the limiting 
case of vanishingly small buckle depths and, further- 
more, that this limiting case establishes a lower bound 
on the critical value of \ for buckled panels. The 
latter conjecture is based on the demonstrated inde- 
on buckle depth for the two-dimensional 


and the more rapid increase of stiffness 


pendence of \ 
problem 

with buckle depth for a finite aspect ratio plate than for 
an infinite aspect ratio plate (beam A test of the 
validity of the conjecture, and others implied in this 
\ppendix, must await a full-fledged dynamic analysis 
ol the aeroelastic behavior of three-dimensional buckled 


panels including large-deflection structural effects 
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The Aerodynamic Action of Triangular 
Horn-Balanced Control Surfaces on the 
Supersonic Delta Wing 


DEREK NAYLOR* 
A. V. Roe & Co., Limited 


SUMMARY 


The pressure distribution of a supersonic delta wing in the re 
gion of influence of a deflected flap type of control with a tri 
angular horn balance is determined on the basis of linearized 
theory. Both subsonic and supersonic leading edges are con 
sidered. General formulas for the control hinge moment incre 
ments AH, AH, are calculated to estimate the effectiveness of 
an aerodynamic balance of the type considered, and correspond 
ing formulas for the increment in the lift AL, are also developed 

The final expressions for the increments are expressed as func 
tions of the delta wing parameter ” = tan IT tan uw (T being the 
leading-edge sweepback and yw the free-stream Mach angle) and 
a geometrical parameter determined by the control surface con 
figuration. These functions are represented graphically. The 
results are applicable to full span flaps, provided the Mach waves 
from the apex of each balance do not intersect the balance area 
or the leading edge of the opposite wing, and to partial span fleps 
if the inboard edge does not lie within the Mach cone from the 
apex of the balance. Formulas for the total hinge moments are 


given in the case of full span flaps 


SYMBOLS 


= wing incidence 


a 

7 = control deflection 

ry = leading-edge sweepback 

m = free-stream Mach angle 

M = free-stream Mach Number 

b = cot yu 

K = tan’ 

n = K/B 

( = delta wing chord (Fig. 1) 

s = flap area behind the hinge line 

b = flap total span (Fig. 1) 

é = mean flap chord (S = bé) 

xo = the maximum forward extension of the balance 
(the length RT in Fig. 1) 

m = Xo/¢ 

x) = basic flap chord (Fig. 1) 

l = Xx)/¢ 

o = m/t 

nN = (1 t m)/(1 t) 

B = (1 +”)/(1 —n) 

o = perturbation potential 

y = pressure function (Eq. 3) 

(26, = rectangular coordinates in the wing plane, with 
x-axis streamwise 

0 = Ky/x 

U = free-stream velocity 

E = complete elliptic integral of second kind with 
modulus n~(n? — 1)1/? 

= = control hinge moment due to wing incidence 

Hy = control hinge moment due to control deflection 


Received July 10, 1956. 
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Now, with the Department of 


AHg, SH, = control hinge moment increments, defined 
Eqs. (6) and (17 

g(o, n = hinge moment increment function defined 
Eq. (9) 

h(r, n = hinge moment increment function defined b 
Eq. (20 

Bc = control lift due to control deflection 

f(a, n = control lift increment function defined by Eq 
(29 


PRESSURE DISTRIBUTION DUE TO CONTROI 
DEFLECTION 


bw TYPE of horn balance considered occupies the 
region RVT of Fig. 1, and it is assumed that the 
hinge line coincides with the flap boundary TT’. To 
find the pressure distribution due to control deflection 
it is sufficient to consider a uniform incidence 7 over 
the control only and zero incidence over the remainder 


of the wing. 


Supersonic Leading Edge 


When the leading edge is supersonic (see Fig. 2 


systematic application of the integral formula,' Eq 


(1), is necessary: 




















Fic. 1 Diagram showing wing geometry 
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denotes the downwash in the wing plane and 
the 


where % 


D y) the interior domain contained between 


Mach lines x — € = +B(y — n) drawn forward from 
the field point (x, y It is possible to apply Eq. (1 


lirectly when the leading edge is supersonic since the 


downwash is known over the entire domain of integra- 
tion. On the control, w = 7, while on the wing, w 0; 
therefore, D(x, y) 
PORS in which the wing region RSS’T may be omitted. 
The configuration of Mach lines shown in Fig. 2 divides 


reduces to the interior of the Figure 


the tip region into zones where the potential takes dif- 
ferent analytic forms. When axes through the tip R 
of the balance are selected, the results for the pressure 
distribution in the regions 1, 2, 3, 4, and 5 influenced 


by the balance are given by 


Wi = (7/2) (1-4 - 
y ys = (1 — n?) cot! [B(x — By) (x + By)]°° 
‘ >) 
¥y;=(1—an '” cot! [B(a By) /(x + By)]'* + (1 2) cos By (x \ 
% = (4/2) + (1 — 2 '”) cot [B(x — By) /(x + By | 7 
where Y(x, vy) = 7Bd,/2Un (3 
; n)/( , ee Se ae eee not known in the region between the leading edge 
, = : ) _ ), ANG Xe GAenovles = axl . or . 
and the Mach wave from the apex R. The influence 


forward extension of the balance (the dimension RT 
of Fig. 1). The overall pressure distribution accord- 
ing to Eqs. (2) is finite and continuous everywhere ex- 
cept across the hinge line, where a discontinuity pl*n B 


in the pressure persists. 


Subsonic Leading Edge 


Eq. 


When the leading edge is subsonic (see Fig. 3), 
1) cannot be applied directly since the downwash is 


y 1 (n + 1)] [(v + By)/(x — Ky)]'* + (1 2) cos 
y lL (n + 1)) } [n(x + By) — (n + 1)xe] (x9 — Ky)} 
y 1 (n + 1)] [@w + By) / (x — Ky)]'* 

vs = (9/2) + [1/(n + 1)] [(e + By)/ (wx — Ky)]' * 


It is noted that TA, one of the Mach lines from the 


hinge axis extremity T (Fig. 3), may meet the sub- 


sonic leading edge and be reflected back across the 


control to determine the zone 2. The condition for 
this is 

nN l+o l 
If, however, l+o2nZ 1 


then the Mach lines from the intersection T intersect 


the supersonic trailing edge and the reflected waves are 


ol no conse¢ juence. 


CALCULATION OF AH 
The hinge moment will be given by 
—2pU Sf (x — x0)o-(x, y) dxdy 
that is, in terms of the function (x, y), 
—4pl? Sf (x — x(x, y) dxdy (5 


evaluated over the control only. 


nll, 


*BH, = 


The hinge moment increment A//, is defined by 


H, — H, (6) 


AH, = 


where /7, denotes the hinge moment exerted on the 


1/2 


io 


[By 


at P of the region RSS’ must then be determined by 
subsonic leading-edge wing tip theory.** In the case 
of a wing of zero thickness the correction simplifies 
Here the method developed by Evvard* is that adopted. 
Evvard's result states that the contributions of the 
zones RSS’ and RR’S’ to the velocity at P are equal 
and opposite so that the integral (1) may be evaluated 
over the region PQR’S’ only. When expressed in 
terms of the pressure function ¥(x, y), the results of the 
integration in the different regions are 


irc =— £ )] 


+ tan—') (x — Ky), [n(: 


control when no horn is present. The hinge moment 
HT, and the hinge moment increment A//, are most con- 
veniently expressed in terms of the functions go, g(¢, ”) 


by writing 
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Mach line system for control with supersonic leading 
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3rk TT, 2pU7c*t® go (7) Sup 
: 7 | 
37K? H, = 2plc*t*® |g) + g(o, n) (S J 
7 1 
so that, from Eq. (6), i twee 
‘ ; 4 
3nK°AH, = 2pU°*c*t® g(o, n (9) / 
— — S7 
To illustrate the derivation of the expressions for kX s' 2 
go, g(a, m) it will be sufficient to set out the steps in the e wr * 
case when the leading edge is supersonic. The results 7 Fy \ ‘\ 
will simply be stated when the leading edge is subsonic. “f J) S \ 
"i rT y, % 
Supersonic Leading Edge ¥ F y, . 
: we i ia oe \Y5|% whe 
When the functions ¥ given by Eq. (2) are inserted y D4 
in Eq. (5) and the integral evaluated over one control, R’ V rs V \ > \ 
the result is F An. \ ; 
B*H./2pUc* (m?/12n) X F, 2 he yala 
}1 + [2 (cos-! n)/ (1 n?)' =}t — [(3 — )t2/6n] (10 / \ \ - 
/ t 
When no horn balance occurs, m7 = 0, so that Q ™ \ . 
_ \ Me ul 
B21, /2pU°c* -[(3 — p)t/6n] (11) \ \ 20 
\ 
Therefore, from Eqs. (7) and (11), a 
i , mi 
go = — [nx(3 — t)/2t] (12 \ lt 0 
Similarly, from Eqs. (6), (10), and (11), Fic. 3. Mach line system for control with subsonic leading edg \ 


B*AH /2ptU/*c? = (m*/12n) X 


Eqs. (12) and (13) apply only when the leading edge is 


,1 + [2 (cos—! n)/r(1 n?)"*]f : Sub 
supersonic, so that 0 TT] | 
so that g(o, v), as defined by Eq. (9), is given in this It 
case by Subsonic Leading Edge 
: 3/9) . z : F 
g(a, n) (na"/2) X In this case, the pressure distribution is given by th cae 
‘ ‘ / . 1/21? 5 . - 5 ‘ ‘ 
» (9/2) + [(cos—! n)/(1 a) "ti §=6( 13 Eqs. (4), and the results of the integrations may bé 
expressed in the form 
29 = |[—nx(3 — nt))/21 for 1 nN 
| 1/2 1.1/2 Tr 
g(o, n) = [no*/(n + 1)|)1 — [(v + 1)/n'*] tan-'n'*} 
(nw/2) |n — 3 — 386 + (2 + 3¢ — o°)n si for | nN | 
f 1/2 1,,1/2) | 
g(o, n) = [no*/(n + 1)] ,1 — [mu +.1)/n'’"] tan-'n'’*} + [n(n — 14+ 6 n ix 
1/2 1/2 9 
[o(n — 1 — o)] —-n “(2+ 380 — o*) tan—'[na/(n | g 
u(t — 3 jo) tan a(n—1—a)] . for | IlI+oSn 16 sucl 
CALCULATION OF A/T, . 
The hinge moment increment A/7/, is defined by equivalent to 
AH, =H, — H 17) ene ae sila Seca : 
h mt*/2E) (3 t l1—t l nN Ze 


i, denoting the hinge moment when no horn occurs, 
Analogously to Eqs. (7) and (8), it is found convenient i ; . 
is : ‘ el ‘ ; lo evaluate control hinge moments with respect to 

to write the moments in terms of functions /rp, h(XA, 2), . mas is elated : 
“oll incidence it is sufficient to consider the pressure dis 
as Iollows: ° : as : ; 
tribution due to uniform incidence a@ over the entire 


3K°H, = pU%c(1 — th (18) wing and integrate the hinge moment of the pressure 
over the relevant control area. To obtain the incré 
2h? - 7273 - 3 Es )] Q) “ ; 
BK* I, pl*c%(I t)® [ho + h(X, n)] (19 ment A//,, it is only necessary to integrate over the 
aun —— | itself. Therefore, 
ch? AH. = pU*ce* — t)* RA, 2) (20) sorn stat Pheretore 
aAH, = —2pU Sf (x — x) o(x, y) dxdy (23 


The formulas for //, are to be found in references 5 and 
6 and, when expressed in the above notation, are integrated over the horn only. 
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TRIANGULAR HORN-BALANCED CONTROL SURFAC 


ove 





Supersonic Leading Edge 
[The Mach wave from the delta apex (see Fig. 4 
livides the wing into zones | and 2. In zone 1, be- 


tween the leading edge and the Mach wave, 
Bo, = Ual l ~ ue 
yone 2, behind the Mach wave, 


Bo Ua(l — n? el (2/8) X 


cos (n- — 6- | f- |! 


when axes through the delta apex are now used and 6 
When the integral of Eq. (25) is evaluated, there are 
two cases to distinguish, depending on whether the 
balance 1s intersected by the Mach wave from the delta 
ipex. 

lf 0 nZjXA <1, the Mach wave does not meet 
the balance and the function /(\, n) defined by Eq. 


20) takes the form 





hi\,n) = n(l—A*A—n 24 i 
; - : Fic. 5. Control hinge moment function o, m) for different 
fO<A <n <1, the wave intersects the balance and values of » 
1,2 12 
\? mmi-—n) “|;l—n [cos~! (Am) + (A8/n) cosh-!(n X)] 
(2 2) [ _ ) = 11/2 ay 2 9- 
(3 + d2) tan n*) /(n® — X + (1 + 3A2) tan! [(1 n*)d* (n \ ; (3 


Subsonic Leading Edge 


In this case, * ‘ 
span, then, in the present notation 


@ la KE) (1 — 6 - . 
S tc-(2 [ a 2 
ind one obtains 
) ( kK, ( tc(2 [ 2 
\,n 1/F) [cos \ + * cosh 1/A) - , 
on ™ . All, L/2)pl"ScACy 
YPA(1 — At) 26 , 
All 1/2 pl ScAC; 


therefore, from Eqs. 20) and (9 


EXPRESSIONS FOR ACy , ACy. : 
7 AC H SC] l 3A (2 [ A(A,n 


Defining the hinge moment coefficients in terms of the 








, : ° s/ Qa ») j i 
mtrol area S behind the hinge line and a mean chord ¢ ACh, Ot 3rA(2 fio, 1 
such that bc where / denotes half the total wing 
Lirt DUE TO CONTROL DEFLECTION 
In this case, it is necessary to evaluate the simpler in 
tegral 
aise nL, Pol SS @Q 
wa 
—T ‘ 
{ aa over the area of the control and wing contained in the 
al - . 
a Mach wave from the apex of the horn balance Phere 
ail . in addition to the pressure distribut the c 
aa 2 fore, in addition to the pressure distribution on the con 
trol, which is given by Eqs. (2) when the leading edge 
is supersonic and by Eqs. (4) when subsonic, it 1s 


necessary to calculate the pressure distribution over 
the region of the wing influenced by the balance. This 


region is zone 5 of Figs. 2 and 3. 
Che results for the pressure distribution in zone 9 





when the leading edge is supersonic may be expressed 








L/e 1} > I? 
a ¥,=(l—n cot~') 3 By + By)}j 
Scheme for finding Ta when the leading edge is super 
SUPe}rsorlic 
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when axes through the apex of the balance are used. 
If the leading edge is subsonic, then 
vs = [1 (n + 1)] [wv + By) (x — Ky)]'? for 1 <n 


Defining the increment for lift in the usual manner, 


AL, = L, — [, 


{(o, n) = nxo*/2 for OS an S 1 
‘ F i > 1/2) 
f(a, n) = (nw/4) [n — 3 — 40 + 2(1 + o)'n } 
i /2 9 ' 1/2 | 
f(o,n) =n’ “(1 + oc)? tan-'[no/(n —1—0)] “+ (n 


USE OF THE CHARTS FOR g(a, 7), h(A, n), f(a, 1) 


Eqs. (13), 
(25), and (26) for A(A, 
calculated, and the results for selected values of » are 
To obtain 


(15), and (16) for g(o, m) and Eqs. (24), 
n) enable these functions to be 


represented graphicaliy in Figs. 5 and 6. 
the hinge moment increments A//,, A/7, for a given 
control and horn balance it is only necessary to com- 


pute the parameters 


g=mit, ’;=(1—t—m)/(1 —2), 2 = K/B 


read off the values of the functions g(o, ”), h(A, m) from 
the charts and substitute in Eqs. (9) and (20). The 
quantities A//,, A//, so obtained then give the incre- 
ments in the hinge moments exerted on a control due 
to an added horn balance of the type considered. If 
the total hinge moments are required, then the func 
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where | a denotes the lift exerted due to control defle 
tion on the control when the balance is absent, and ] 
denotes the lift exerted on the control and the win, 
regions influenced by the balance, and introducing 


function f(o, ”) by means of 
wh °AL, = 2pt/*c*t* fle, 2) 2 


one finds 


for l1iouQgl+ea 
5) 
. » 1/2 
2) (n — 3 — 4c) tan-' [o6 (n — 1 — o)] ~ - 
‘ 1/2 ° 
(n/2) [o(n — 1 — o)] for | l+oa n 


tions go, 4y may be calculated from Eqs. (12), (14 
(21), or (22) and Eqs. (8) and (19) may be used, pr 
vided full span flaps are being considered. In th 
case of partial span flaps, the results are still applicabk 
as here derived provided the inboard edge does not li 
within the Mach cone from the apex of the near balance 
and, in this case, the quantities A//,, A//, may b 
added onto data for such flaps without the hon 
balances. 

Similarly, Fig. 7 gives the values of the functio 
f(o, n) computed from Eqs. (30). AZ, may then | 
obtained from Eq. (29). Reference should be mad 
to the work of Coale,’ 
horn balance mounted on 
parison of the results is difficult, as both the type \ 
balance and the basic flap differ from those consider 
The geometrical configuratior 


rectangular 
Con 


who considers a 
a rectangular flap. 


in the present work. 


Continued on page 610 
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On 


roducing 


‘ixing Boundary-Laver Transition 
Supersonic Wind-Tunnel Models 


2) MARVIN LUTHER* 


California Institute of Technology 


SUMMARY 


to reduce or eliminate scale effects of supersonic W1 id 
7 dels, it often is important to fix transition to turbulent 
the boundary laver of the model, by some artificial 
\ systematic series of wind-tunnel tests was conducte 
1? 14 id e use of single-element and distributed roughness for 
urpose, and to determine the character of roughnesses 
used, pri in SCS a 
! I would promote transition with minimum collateral effeet on 
, n th flow. Curves are presented showing critical roughness height 
‘ qnre ] 
applicabl function of Mach and Reynolds Numbers; also presented are 
oes not li schlieren photographs showing the effect of Reynolds Number, 
ir balance roughne height, and angle of attack on the transition point 
: models used for the study are a family of cone-ogive-cvlinders 
a may b : 
; fineness ratio of 10.6 
the hor 
. F SyMBOLS 
functor 
v then b linder cross-sectional area vase area 2.011 
be mad ci 
model lateral surface area 


‘ctangular ape 
2 model lateral surface area forward from station ¢ 


Pp. Com constant defined by Eq 6 
le type ( ( lrag coethicient 
‘-onsideret ( drag coefficient for model with laminar boundary laver 
figuratior drag coetiicient for model with turbulent boundary 
i\ 
drag foree measured by balance 
eight of roughness element 
ee | Mach Number 


tatic pressure 
dvnam 


Revnolds Number 


pressure pur/ Zz 


Sutherland’s constant = 216°R 

static temperature 

velocit 

oordinate along model surface, positive aft, with ori 
gin at tip of model 
rdinate normal to model surface with origin at 
model surface 


indicated angle of attack 

ti of specific heats 
boundary-layer displacement thickness 
VISCOSITN 
kinematic viscosity 


iss density 





Subscripts 
adiabatic surface condition 
1 odel b is¢ 
critical 
= station at which the roughness element is located 
transition point 
odel surface 


free-stream conditions 


ee Received December 26, 1956 

d Formerly, Research Engineer, Jet Propulsion Laboratory 
Now, Research Engineer, Military Operations Research Division, 
Lockheed Aircraft Corporation, Burbank 
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1) INTRODUCTION 


S* supersonic wind tunnels are limited in size and 


stagnation pressure by power requirements and 


other considerations, many model tests are made at 


Reynolds Numbers much lower than those of the full 
As a the 


layer may be laminar over a large part or all of the 


scale aircraft or missile result boundary 
model surface instead of turbulent as it frequently is at 
the higher Reynolds Numbers of flight. This difference 
the rate of growth of the 
the 


the conditions under which the boundary 


affects the skin-friction drag 


therefore distribu 


and pressure 


boundary layer 
tion lave 
separates, its behavior after separation and the effec 
tiveness of control surfaces that are partially immersed 
in the boundary layer. As the flight speed increases 
these effects become more important because of the fact 
that, in general, the boundary-layer thickness increases 
with Mach Number (at constant Reynolds Number 

In order to reduce or eliminate these scale effects it is 
often necessary to fix transition to turbulent flow, 1n the 
boundary layer on the model, by some artificial method 
This requirement presents the problem of choosing a 
mechanism which will promote transition with minimum 
collateral effect on the flow 

Although a fair amount of data on this subject has 
been published, cf. also, bibliography of reference | 
most of the material concerns two-dimensional models 
and or low-speed flow. A series of five tests were pe! 
formed! at the Jet Propulsion Laboratory during the 
November to 


period from December, 1954 to 1955 
study the production of transition on three-dimensional 
models at supersonic speeds; the tests were conducted 
in the 20-in. wind tunnel on cone-ogive-cvlinder models 
at Mach Numbers 1.63, 3.07, 4.09 ».04 The 
method used to produce transition was the addition of 


and 


roughness to the tip of the cones, two types of roughness 


being used—-single-element (wire) and distributed (car 


borundum grit 


2) NMIETHODS OF FIXING TRANSITION 


The methods which have been used for fixing transi 


tion? include (1) attaching carborundum grit, single or 


multiple wires, or ridges to the model surface; (2) cut 


ting multiple grooves into the surface; and (3) injecting 


air into the boundary laver. The attachment of grit or 


+ The term “‘fixing transition” is used to indicate the artificial 


establishment of the boundary-layer transition point at or near 


the device used for that purpost 
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windshield 


Jasic model configurations with 


Fic. | 


wires to the model surface appears to be the most prac 
tical technique for use in supersonic wind tunnels, and 
the present investigation was limited to the use of grit 
and single wires. Multiple wires may prove more satis 
factory in some situations since their size can be smaller 
than that of the single wire, and consequently the dis 
turbances produced in the free stream might be ex 
pected to be less severe. Inasmuch as grooves have 
been found less satisfactory than ridges,*’ they were 
not considered in this study. Air injection is probably 
the most elegant way of producing transition, and the 
size of the disturbance can be varied by controlling the 
air flow. Its application to general models is limited, 
however, because of the increased complexity of the 
models. 

Carborundum grit has the disadvantage that it tends 
to be eroded away in time. If grit is used in a lengthy 
test, the forward tip of the model must be readily de- 


CAL SCIENCES—AUGUST, 1957 
FABLE | 
Roughness Elements 
Mean Grit \merican \\ 
Grit No Height (in Wire Gage No. Diamet 
320 QO. OO14 3S 
2() 0. OOL2T 333 
150 0) OOB9 4 ) ‘ 
120 0.0052 2Q 
100 0 0066 a7 ) 
S() Q OO94 24 
Hu O O18] 14 ) 


tachable in order to permit reapplication t the grit 


the end of the day. A possible advantage of the grit ; 


that, in low-speed tests, it has been shown to provi 


elocit 


less distortion of the turbulent boundary-laye1 


profile. The average size of grit used during the test 
are listed in Table 1. The sizes were supplied by th 
manufacturer and correlate well with measurements 
made with an optical microscope. The sizes of wir 
used as single element roughnesses are also given j 
lable 1. 

(3) DESCRIPTION OF MODELS 


The models were bodies of revolution composed 


conical noses and cylindrical afterbodies joined by tar 


gent circular ogives. The three basic models had small 
radius, tangent, circular ogives and cones of 12, 20, an 
30 degrees included angle; they are designated models 


NB), NooB;, and NzoB, 


uses the same conical tip as NB, and the same after 


respectively. A fourth mod 


body, but has a large-radius ogive his model was 
designed for the purpose of observing the effect of pres 
iB Che four 


dimensions Fig 


sure gradient and is designated N 


models are shown with basic iN 


Model profile coordinates are given in Table 2 
Several tips were made for each cone in order to facili 
tate the exchange of roughness elements. Nose tips {01 


the 20° cone, with typical roughness elements applie 


» 


Coordinates 





TABLI 
Model Profile 
N; N 
Station fT Diameter Station Diameter 
(in. ) in in in 

0.000 0.0000 0.000 0.0000 
7.3809 1.5351 +. 160 1.4693 
7.399 1.5450 +.210 1.4861 
7.409 1.5531 $+. 260 1. 5020 
7.459 1.5611 $.310 1.5170 
7.509 1.5677 $ 360 1.5307 
7.559 1.573 $.410 1.5430 
7.609 1.5787 $. 460 1.5542 
7.659 1. 5835 $.510 1. 5641 
7.709 1. 5876 $. 560 1.5728 
7.759 1.5917 $610 1. 5800 
7.809 1.5948 $ 660 1. 5860 
7.859 1.5975 t.710 1.5905 
7.909 1.5991 +. 760 1. 5941 
7.959 1. 5996 $810 1.5965 
8. 459 1. 6000 $ 860 1. 5982 
& 959 1 6000 $910 1.5991 

$. 960 1.5998 

5.460 1. 6000 

§ 950 1.6000 


+ Station numbers in this report 


are given in inches af 


Noo N 

Station Diameter Station Diameter 

in inl in in 
0.000 0.0000 0.000 QO OOO 
2? OOO 0.7021 * Seo 1. 4468 
2.250 0.7827 2 762 1 4677 
2 SOO 0.8590 » £12 1.4868 
2.750 0.9323 » S62 1. 5040 
3.000 1.0029 2.912 l S995 
3.250 1.0695 » 962 1.53 
3.500 1.1333 3.012 1.5468 
3.750 1.1942 3.062 1.5582 
£000 1.2491 $.1I2 1. 567¢ 

$250 1.2978 3.162 1. 575¢ 

$ 500 1.34235 3. 212 1.5822 
$4. 750 1. 3854 3.262 1.5878 
5.000 1.4252 3.312 lL. 092 
5.250 1.4602 3.362 1.596 
5.500 1.4900 3.412 1 5987 
5.750 1.5150 3.462 1. 600: 
6. O00 1. 5865 & 962 1. 6006 
6.500 1.5718 
7.000 1. 5889 
7.700 1.5986 
8.960 1. 6000 


from the theoretical tip before the tip radius was applied 
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g the tests 


re shown in Fig. 2. The tips were slightly blunted, 
with tip radi ranging from 0.0025 to 0.0067 body 
jameters. It is believed that this degree of bluntness 
id not have a measurable effect on transition 

The base of the models was reamed out for a depth of 
, in. so as to leave a 0.081-in. wall thickness. The 
windshield was provided with a collar of such diameter 
sto leave a 0.051-in. annular opening between wind 
nd model base. The reasons for this design and 
etails of the base-pressure instrumentation are dis 
ussed in Section 5 


!) WInpb TUNNEL AND INSTRUMENTATION 


Design features of the 20-in. wind tunnel permit its 
tinuous operation in either a variable-density closed 
ireuit or an open circuit. For the present tests it was 
perated as a closed-circuit tunnel. Any value of test 
section Mach Number between 1.33 and 5.04 may be 
set by remote control from the control-room console 
Normal operation is restricted to fourteen Mach Num 
bers, between 1.53 and 4.54, for which flow surveys have 
wen made; however, three runs were made during the 
earlier tests at Mach Number 5.04. It is expected that 
idditional runs will be made at 5.04 when planned tun 
nel modification and calibration make operation at this 
Mach Number more convenient and reliable. The tur 
bulence level is approximately 0.005, defined as the 
ratio of the disturbance velocity in the streamwise 
direction to the stream velocity, and measured in the 
supply section. Data from flow surveys for the Mach 
Numbers used are presented in Table 5, together with 
other test-section flow parameters; at all Mach Num- 
bers the test section is 20 in. high by IS in. wide. Dur 
ing a given run the air supply is maintained at a con 
stant temperature and within a dew-point range of —50 

10°F. Test-section flow parameters were deter 
mined by using the Mach Number, the measured sup 


ply-section stagnation pressure and temperature, and 


the assumption of isentropic flow in the nozzle. 

Of the various methods which have been used to de 
termine the location of transition, those which were best 
idapted to the present test were (1) the measurement of 


lrag, and (2) schlieren photography. 








—— 





Model tips with typical roughness elements 


PABLe 3 
Pest-Section Flow Condition 


Potal pressure (psia 


Maximum 24.1 29 9 52.0 65.5 

Minimun 1.9 >.4 r 7.0 
Dynamic pressure (psi 

Maximum | $} NS 

Minimun 8 25 é , 4 
Static pressure (psia 

Maximum 54 

Minimum 4 , 4 O5 
Reynolds No./in. & 10 

Maximum ) 5 13 S l 

Minimun O5 IS 0.1 } 


Mach Number variation FO OO t Olt +0) 


Flow-angle variation +) 2 +1) 


stat essure Variati 
Af | ) >. 
Flow calibration has not bee de at this M r 
Phe models were sting-supported from a crescent 


shaped strut carrier, which may be rotated 30° in pitch 
The sting and strut carrier are protected from aerody 
namic loads by a windshield which is open to the tank 
surrounding the external balance Flow between the 
test section and the balance tank is prevented by evacu 
ating the tank to the model base pressure The base 
pressure used for this purpose was that measured at an 
orifice located at the periphery of the collar on the wind 
shield, since this pressure should be very nearly equal to 
that at the annular orifice between the windshield collar 
and the model base. A second pressure, that in the 
chamber at the model base, was measured on a mercury 
micromanometer which could be read to 0.001 em. of 
mercury. This pressure, which might differ from the 
base pressure in the advent of a small flow through the 
model base, was used in the data reduction 

The balance system mechanically separates the model 
lift, drag, and side forces and the pitching, rolling, 
and yawing moments. These forces and moments are 
sensed by hydraulic load cells and are measured by self 
balancing beams. The moments as measured are re 
ferred to a set of orthogonal axes which intersect in a 
point that is coincident with the center of rotation of 
the support system and the center line of the test-sec 
tion windows. 

Flow in the test section may be observed optically by 
either the shadowgraph or the schlieren technique. For 
the first four tests the schlieren system was used with a 
spark source of approximately %-microsec. duration. 
For the fifth test a l-microsec. source was used and was 
found to be less satisfactory than the source used for the 


earlier tests. 


)) Test PROCEDURI 


During each run, data were obtained at four to six 
stagnation pressures, including the maximum and 
minimum available, in order to provide variation of 
Reynolds Number. At each pressure level the drag 
force and the base pressure were recorded and reduced 


to coefficient form at l-min. intervals until the drag co 
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Fic. 3 ypical plot of drag coefficient vs. Reynolds Number 


efficient repeated within the estimated accuracy of the 
system. Assuming, then, that equilibrium conditions 
had been reached, the pressure was changed to the next 
level. A spark-schlieren photograph was taken at each 
point—-i.e., each time the drag data were recorded. A 
period of 2 to 6 min. was required for the complete pro 
cedure at each pressure. Plots of drag coefficient vs. 
Reynolds Number were made as the run progressed, and 
any data which did not appear reasonable were checked 
before the conclusion of the run. 

It was noted that the data were very sensitive to the 
accuracy with which the base pressure and balance-tank 
pressure were balanced. Accordingly, the pressures 
were balanced within | mm. of mercury when possible. 
At high Mach Number and low supply pressure (or 
Reynolds Number) the base pressure was lower than 
the lowest balance-tank pressure because of limitations 
of the pumping system; when this condition obtained, 
the measured drag data were questionable. The data 
were also aifected by the condition of the model surface; 
as a result, the surface had to be carefully cleaned after 
each model change. 

During several runs the model was pitched, and 
schlieren photographs were made at several angles of 
attack; however, drag data were obtained only at zero 


angle of attack. 


(6) DATA REDUCTION 


The drag data were reduced to coefficent form using 


the relation 
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~ 
~ 
~ 
| 


hydray 


where PD is the drag force recorded by the 
balance, /, is the base pressure (within the base cha; 
ber), A is the base area (equal for these models to t 
reference area A), and fp and q are free-stream stay 
and dynamic pressure, respectively. 

The drag coefficients were plotted vs. Reynolds ny 
ber Rey and then crossplotted vs. roughness height 
at constant Revnolds Numbers. The critical roughn 
height was inferred from these plots; Figs. 3 and 4 
typical. To aid the fairing and intepretation oj 
drag-coefficient plots, the transition-point locatior 
was measured from the schlieren photographs 
plotted vs. Reynolds Number. The same quantity y 


also computed from the drag data using the relatio 


A; A Cp ~ Cpr Co: Cu 


where A, is the area of the surface forward from stati 


v, and Cp, and Cpr are the drag coefficients for tl 


model with a laminar and a turbulent boundary layer 


respectively. The values of Cp; and Cpr were obtain 
by calculating approximate values, plotting these valu 


vs. Reynolds Number, and shifting these curves to mal 


them conform to the values measured for models whic! 


from the schlieren evidence, had completely laminar or 


completely turbulent boundary layers. 


To attempt correlation with Dryden's roughnes 


parameter” (k 6,*), it was necessary to compute th 


boundary-layer displacement thickness 6,* at the st 


tions where roughness was applied to the conical tips 


For this purpose, the defining equation 


” ple Coll iy 
) 


was used, as was also the solution to the flat-plat 
boundary-layer problem given by Chapman and Rube 


sin’ which reduces Eq. (4) to the form 


5" (vexC/uy)’’” | 1.7208 1.1094(4 1). 
1.9384] (t,. ¢ (tar / to) }j 
‘ 
4] 2 
u a ’ | es RES VE Cae eocaail 
=) = cl +—0F Se | 








Fic. 4. Typical plot of drag coefficient vs. roughness 


height 
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where /, is the wall temperature, ¢,,, is the adiabatic wall 
temperature, and w, Mm, J/), and ¢ are free-stream 
kinematic viscosity, velocity, Mach Number, and static 
temperature, respectively. The quantity C is a con- 
stant defined by the assumed relation 


b/g = Cl(t/ to) (6) 


where » and wy are the local and the free-stream vis 
cosity. This relation provides reasonable accuracy at 
moderate Mach Numbers when the constant C is de- 
termined by substitution in Sutherland's formula, 


where S is Sutherland's constant (216°R. for air). 

The cone solution was obtained using the transforma 
tion law which states that the laminar cone solution is 
the flat-plate solution computed for one third of the 
Reynolds Number on the cone.*’ The “‘free-stream” 
quantities in Eqs. (5) and (6) are computed on the cone 
surface, behind the conical shock. 

For an insulated model the last member of the right 
hand term of Eq. (5) becomes zero. This approxima 
tion was made for the models used in the present tests. 
Comparison of the computed displacement thickness 
with the apparent boundary-layer thickness measured 
on the schlieren negatives showed agreement within the 
limit of resolution of the film (approx. 0.002 in. 


(7) RESULTS 


(A) Critical Roughness Height 

The principal purpoce of the tests described herein 
was to obtain the critical roughness-element height &, 
(defined as the minimum height sufficient to produce 
transition at or very near the roughness element) as a 
function of model and flow parameters. Figs. 5 to S 
present plots of k, vs. Mach Number and Reynolds 
Number per inch for distributed and single-element 
roughnesses. In these Figures the dashed portions of 
curves indicate either extrapolations from the measured 
data or curves based on too few or uncertain points; 
they should be used with this fact in mind. 

A comparison of these figures shows a rapid increase of 
k. with increase of Mach Number, particularly at low 
Reynolds Numbers, as would be expected because of the 
increase of boundary-layer thickness. It is wso seen 
that, for the same roughness height at the saine station 
x,, the distributed roughness is more effective than the 
single-element roughness, and that the margin o! eifec 
tiveness decreases with increase of Mach Number. The 
difference in effectiveness is not surprising since the 
distributed roughness extends forward from station x 
to the tip of the cone. 

The effect of model geometry can also be seen. At 
Mach Number 4.09, both the data for distributed rough 
ness and for single-element roughness indicate that k 
tends to increase with increase of cone angle, especially 
at the lower Reynolds Numbers, and that the effect of 
ogival radius on k, is negligible. At Mach Number 1.65 
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neither of these trends is consistently observed: hoy 
ever, the spread of the curves is relatively small 

To check the effect of the amount of distribu 
roughness used, data were obtained for model .V.B 
Mach Number 4.09 using grit number 320 with 4 
O.1S75 and 0.500 in. Within the accuracy of th 
measurements, it was found that covering the add 
tional surface area with roughness did not further a} 
the character of the boundary layer. 


(B) Transition Point 


The transition point was obtained from spark 
schlieren photographs and from the drag data [cf. Se 
tion (6)]. On the average, the drag data indicated 
slightly earlier transition than did the optical data, sug 
gesting that the point observed on the schlieren phot 
graphs was aft from the transition region (in which th 
shear stress at the model surface is higher than in th 
laminar region). 

Typical spark-schlieren photographs showing th 
transition point are presented in Figs. 9 to 12. Th 
effect of Reynolds Number on transition is illustrate 
Fig. 9, and the effect of grit height is shown in Fig. || 


As the grit height is increased, the transition point 


moves forward to a point near the ogive and remains 
fairly constant as the grit height is increased further 


until a near-critical value is reached, at which time the 
transition point moves onto the nose. It is believe 
that, in the absence of the negative pressure gradient 
over the ogive, the transition point for the models with 
120 and SO grit would be farther forward; it 1s further 
believed that the influence of the pressure gradient 
increases as the transition point approaches the ogive 

The effect of angle of attack on models with an 
without grit is shown in Figs. 11 and 12. As would be 
expected, the transition pont on the lower surfa¢ 
moves aft and that on the upper suriace moves forwari 


with angle of attack, these movements being much less 


pronounced on the model with roughness added. 

To correlate the data on fixing transition by the addi 
tion of roughness elements, Dryden” suggested the ust 
of the parameter & 6,*, where & is the roughness height 
and 6,* is the boundary-layer displacement thickness al 
the location of the roughness element; this parameter 
provided excellent correlation of data for low-speed flow 

\J = 0) with zero pressure gradient. It was found bi 
Goddard? and by Brinich* that data for supersoni 
Mach Numbers did not correlate well with the rough 


ness parameter, yielding curves that were of similat 


shape but with values of & 6,* several times larger 


(Goddard's data indicate that to obtain a single-curvt 


correlation would require modification of the roughness 


parameter by additional functions of Mach Number 
Reynolds Number, and roughness size. 

Plots of k 6,” 
Numbers with and without roughness were prepare 
using data from the present tests; they are not pr 

j Unpublished data obtained by Frank E. Goddard, Jr., at! 
Jet Propulsion Laboratory, C.I.T., on cones with single-elen 


roughmnesses 


vs. the ratio of the transition Reynold 
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sented in this report, however, since the roughness 

parameter failed to correlate the data. The curves had 
the same character as those obtained by Goddard, with 
the pressure gradient over the ogive as an additional 


variable. 


(S) CONCLUSION 


Data have been presented showing the effect on 
boundary-layer transition of adding single-element or 
distributed roughnesses near the tips of cone-ogive 
evlinder models, using roughness height, cone angle, 
ogival radius, Mach Number, and Reynolds Number 
Each model had a fineness ratio of 10.6. 
Mach 
and Reynolds Numbers on the critical roughness height, 


as variables. 
Curves are included which show the effect of 
defined as the minimum height sufficient to fix transi 
tion at the element. 

It is concluded that, for the ranges of variables tested, 
both 


nesses can be used satisfactorily to fix transition, with 


the single-element and the distributed rough 
negligible efect due to the roughness itself on drag 
(apart from that attributed to movement of the transi 
that 
roughness size may have such an effect on the wave 


drag, especially at lower Mach Numbers. 


tion point It is believed further increase in 
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949 ved in dequate Recent large deflection theoretical treat 
undary Lay ents of the problem have shown reasonable correlation with Statement of Problem 
gle I ttack periments but require prior knowledge of the initial imper 
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inalyzed semiempirically In the analysis best-fit curves 
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re presented using theoretical parameters and shapes of curves ; : ite. 
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mpressive buckling curves are then developed as 90 per cent structures. In addition to obvious manufacturing ad 
bability curves from the test data. In general, these statis vantages, the monocoque type of construction reduces 
tically defir sig ves sienific » lowe : ° 7 i 
Ry dcined Genge carves are sigmmcantiy lower than prev thermal gradients and is very well suited for integral 
islv available design curves . : mat “ 
fuel tanks. Ring-stiffened shells are a type of con 
ena struction intermediate between fully stiffened and fully 
. 30LS . 
monocoque. If the rings are sufficiently rigid to pre 
— et ee ee Cee. & vent general instability, the portions of the shell be 
CEtt ; ike. 
cating Th } viii ll atin tween rings may be considered as monocoque with ap 
Tho Ss of Clastic ’ ° 
mic Spe K buckling stress coeflicient in the equation «a... 7 propriate end conditions. rhe present investigation 


K,|x?E/12(1 — y?)](t?/L? 

wide plate column buckling coefficient 

length of evlinder 

internal pressure 

rif radius of cylinder 

Institute wall thickness of evlinder 

val Resear , > 
/ curvature parameter = (L? 

plasticity correction factor 


Poisson’s ratio (assumed equal to 0.3 


o unpressurized critical buckling stress 
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* Formerly Structures Engineer, Missile Develop 


nt Now, Sei 
tories, Cook Electric Company 


Senior 


ior Aeronautical Engineer, Cook Research Lab 


was initiated to develop reliable design information on 
the critical axial compressive buckling stress for thin 
walled circular cylindrical shells with and without in 
which considers axial 


ternal This 


loading only, is the first part of a general study covering 


pressure. paper, 
the buckling of circular cylinders and curved panels 


under compression, bending, torsion, and combined 


loads. 

The buckling of many structural elements occurs in a 
smooth wave form which often corresponds closely to 
classical theory 


the sinusoidal 


Cylinders, on the other hand, often snap-buckle into a 


shapes assumed in 


sharply crested large-deflection pattern similar to a 
In this large 


diamond pattern can be seen to grow out of some initial 


series of diamonds. some instances 
irregularity of the curved surface. 

It is well known that for unpressurized evlinders large 
discrepancies exist between values of the critical buck 
ling stress determined by test and those predicted by 
theory. 


small-deflection Experimental 


the classical, 
values of the critical stress for unpressurized cylinders 


are always lower than those calculated theoretically 
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and are in some instances only about 10 per cent of the 
predicted values. Large-deflection theories have been 
developed by von Karman and Tsien' and others* 

in an attempt to account for this discrepancy. These 
large-deflection theories explain the so-called snap 
buckling phenomena by showing the existence of a large 
deflection equilibrium configuration at a lower stress 
than that computed by the small-deflection theory. 
Nevertheless, the critical stresses computed by these 


theories are greater than those found experimentally. 


Another explanation of the behavior of unpressurized 
cylinders in compression has been the derivation by 
Donnell and Wan? of a large-deflection theory which in 
cludes the effects of initial irregularities and accounts 
for the plastic behavior of the shell. This latter ap 
proach appears to offer the best correlation with 
experimental results where the initial eccentricities 
were known. There results also suggest the depend 
ence of the critical stress coeflicient on the 7 ¢ ratio be 
cause the initial imperfections of an unstiffened cylinder 
are likely to be greater as the 7 ¢ ratio increases. In 
spite of the theoretical soundness of a large-deflection 
theory which includes the effects of initial imperfec 
tions, the procedure has the disadvantage that it re 
quires a prior knowledge of the initial imperfection. 

A discrepancy also exists between theory and test for 
pressurized cylinders. The compressive buckling of 
pressurized cylinders was investigated by Fliigge? who 
used small-deflection theory. His conclusion that the 
buckling stress was not aifected by internal pressure is 
in contrad.ction to experimental results. Based on von 
Karman and Tsien’s theory, Lo, Crate, and Schwartz' 
have derived a large-deflection analysis for pressurized 
cylinders. This large-deflection theory indicates an 
increase in the buckling stress with pressure until a 
maximum value is reached, but the experimental values 
are lower than those predicted by theory. This dis 
crepancy, although not as large as for unpressurized 


cylinders, is likewise thought to be a result of the init 


imperfections of the test specimens. 


Scope of Investigation 


As pointed out above, theory does not adequately 
scribe the compressive buckling behavior of circul 
shells without a prior knowledge of an initial imperfi 
tion parameter. In the authors’ opinion, it would | 


difficult, if possible, for designers or analysts to deter 


mine an initial imperfection factor for aircraft stry 
tures. Therefore, it has been necessary to resort t 


semiempirical development of design curves. Th 


semiempirical method has the disadvantage that ther 


is considerable scatter in the experimental data. Hoy 
ever, this approach is thought to offer the most practi 
solution to the problem at the present time. The bas; 
assumption in this semiempirical method 1s that th 
initial imperfections of the test cylinders were repr 
sentative of those likely to be encountered in actu 


aircraft structures. 


The theoretical parameters derived from Donnell’ 


linear small-deflection theory’ have been used in th 
analysis of the unpressurized cylinders. In addition t 


avoiding a direct prediction of the initial imperfections 


the use of the parameters from Donnell’s equation has 


the advantage of allowing an analysis over the wid 
range of geometry from the flat plate to the long cyli 
der. Once the theoretical parameters had been choser 
the test data were analyzed using the shape of the 
theoretical curve as a guide. A best-fit curve 1s fitte 
to the test data and statistical design curves are derive 
on the basis of pr Ibability theory. 

For pressurized evlinders, parameters from the theor 
of Lo, Crate, and Schwartz® have been used in the di 
termination of an empirical design curve. This pr 
cedure avoids the necessity of predicting the effect 
initial imperfections. A design curve based on prot 
ability theory has been drawn indicating the increas 
in the critical stress as a function of a pressure paran 
eter to facilitate design interpretation. 

In reviewing the prior experimental investigations 
became apparent that a considerable gap existed 1 
the available test data, especially in the tests of pres 
surized cylinders. In order to overcome this lack 
data, an extensive series of tests was performed over 
range of parameters likely to be encountered in aircralt 


structures. 


EXPERIMENTAL INVESTIGATIONS 


North American Aviation Investigation 


An extensive series of compression tests of unpres 
surized and pressurized model circular cylinders was 
performed at the Missile Development Division Stru 
tural Test Laboratory. At the initiation of the tes 
program in 1951 only a small amount of data over lim 
ited ranges of the parameters was available on t! 
compressive buckling of pressurized cylinders.  Becaus 


of the scatter in compression tests of unstiffened cyli! 
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lers, it was deemed necessary to perform an experi 
mental investigation which would extend the test data 
within the range likely to be encountered in missile or 
ircrait structures. 

\ total of S83. tests was performed, 14 of unpressur 
ed evlinders and 69 of pressurized cylinders. The 
test evlinders were fabricated from 0.005-in.-thick 2S 
H-IS aluminum foil and 0.0032-, 0.004-, and 0.00S7-in. 
thick half-hard 1S-S stainless-steel foil. The above 
sheet thicknesses were the averages of a number of 
micrometer measurements. The cylinders were pre 
pared with lengths of 10.5 or 21.5 in. and a radius of 
S79 in. and had f ratios of approximately 1,000 


750, 2.190, and 2,730. The internal pressure was 
Most of the cylinders 


) 


varied from 0 to about 24 psi. 
were fabricated by wrapping the sheet around a form 
ind adhesive bonding a single 3 S-in.-wide longitudinal 
seam with Epon VI adhesive. Some of the later 
evlinders were fabricated by seam welding, which 
proved to be a faster and simpler technique 

A schematic drawing of the test jig with a cylinder 
in position for testing is shown in Fig. 1. The test 
cylinders snugly fitted the lower and upper heads of the 
test jig and were clamped to the heads by metal straps. 
[he straps prevented slippage and leakage between the 
heads and the cylinder. Internal pressure was pneu 
matically supplied, but for reasons of safety the cylin 
ders were partially filled with water for all tests in 
which the internal pressure exceeded 2 psi. A mercury 
manometer was used to measure the internal pressure. 
\n internal pressure equivalent to a head of water ot 
half the specimen length was added to the applied 
neumatic pressure for the analysis of those specimens 
which were filled with water. The axial compression 


loading 


was applied by a hydraulic strut and was 
measured by an SR-4 load cell. 

Phe geometrical parameters and the buckling stresses 
of the test evlinders are recorded in Table 1. The 
buckling loads of the test specimens were determined 
visually. Buckling occurred at or just prior to the 
ultimate load the specimen would resist. For most 
specimens visual determination of the buckling load 
was relatively simple because the buckles often snapped 
into position. Occasionally, buckles appeared to grow 
out of relatively large initial irregularities. 

Photographs of typical buckle patterns in unpres 
surized and pressurized cylinders are shown in Figs. 2 
ind 3. These diamond-shaped buckle patterns are 
typical for compressive buckling in unstiffened cylin 
ders. Internal pressure changed the buckle pattern 
by decreasing the axial and circumferential wavelengths, 
though not in the same proportion. The buckle pat 
terns frequently ran through the vertical seam, but in 
no case did buckling begin at the seam. 

Incipient buckling of a cylinder was often indicated 
by the formation of small circumferential ripples. Al 
though very shallow, these ripples (similar to the shape 
of an axisymmetrical compression buckle) could easily 
be seen in the reflecting surface of the cylinder as shown 


in Fig. 4. Once the specimen was close to the buckling 
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TABLI 
Resuits of North American Aviation Tests 
r psi psi 5 8 

2S-H18 a linders (E x 10° 8.75 
35(1) |21.5 7 855 33 1 1.9 2,12 
36(1) ‘ ‘ 818 38(1) i i 2.6 2,14 
40(1) 1,51 39(1) 2.¢ 2,18 
41(1) 1,782 34(1 ' ' 3 818 
42(1) ’ ’ 1,618 3 21.5 5 - 
32(1) (21.5 | 175 9 2,073 

8-8 1/2 hard s ess steel 1 ers (E 27 x 8 
43(1) /21.5 | 1006 5,81¢ 23(1) 2734! 8 4 
44(1) ‘ ‘ . 6,57 31(1) i i 8.2 4 
45(1) ' 5,80C 28(1) 24.2 3, 
56(1) 6,78 30(1) 24.2 3,2 
303 1.0 8,946 (1) ’ 24 4,2 
304 1.0 8,654 1(1) 9.5 2 4,3 
46(1) 4.4 9,47¢ 69(1) | 21.5 4,063 
47(1) ‘ 93 2(1) i 4.4 4,347 
55(1) ! 6,15 2(2) i 4,1 
57(1) 4 2(3) 3,324 
48(1) 6 3(5) 4,233 
491) 8 9(1) 4,4¢ 
501) 8 66(1) 4,29 
305 68( ’ 4,4¢ 
51(1) é 2(4) 8.4 3,35 
52(2) é 2(5) i 3,864 
53(1) 2 3(¢ 5,114 
54(1) 2( 3(8 4,205 
214(1) ® | 24 6( 4,2 
215(1) é 4 8(1 ’ 4,2 
231(1)) 2188 ' 9(2) 8.4 4 
232(1)' 21.5 2188 24.4 70(1 0.4 4,631 
25(1) 9.5 | 2734 71Ql 4 4,744 
29(1) 9.5 ‘ ‘ 2(6) 6.4 3,92 
4(1) [21.5 7€ 3(9) j 3, 
4(2) ‘ ¢ ’ 4,318 
71 1,097 6.4 4 
7(2) l, 65(1) 20.4 4 
58(1) ' ’ 1 3( 24.4 3,58 
é ) 21.5 0 1 4(3 4 4,034 
19(1) 9.5 4.2 4, 13 3,9 
20(1) ‘ 4.2 4,375 14(1) 3,8 
22(1) 4.2 4,21 15 4 
16(1) 8.2 4,489 64( ’ 4,09 
17(1) ’ 1 8 4,37 2 ’ ' 4.4 3,63 
21(1) 9 2734 . 4 ¢ 34 ‘ 


load, a buckle could be easily induced by jarring or by 
a small lateral finger pressure. A few test specimens 
were discarded because buckling was precipitated by 
finger pressure which caused the buckle pattern to snap 


suddenly into the evlinder 


Summary of Other Investigations 


An important aspect of the semiempirical analysis 
has been the proper accumulation and utilization of test 
data. A large number of investigations were carefully 
reviewed to obtain test data, but only the experimental 
results from the eleven investigations’ ? summarized 
in Table 2 have been used in developing the analysis 
he results of other investigations were not included 
in the study for one or more of the following reasons 
Insuflicient data were given; test specimens were fabri 
cated from nonmetallic materials such as paper, cellu 
loid, or rubber which frequently exhibit variations in 
properties and thickness; critical stresses were above the 
vield point; critical stresses were above the propor 
tional limit and insufficient data were given to make a 
plasticity correction; the fabrication procedure was not 
rauio was 


typical of aircraft construction; and the » 


less than 100. 


ANALYSIS OF UNPRESSURIZED CYLINDER DATA 


Theoretical Considerations 


The compressive buckling behavior of unpressurized 


cylinders may be considered in four classifications 





TOURNAL 


Fic. 2 (top). Typical buckle patterns for unpressurized cylinders 


Fic. 3 (bottom). Typical buckle patterns for pressurized evlinders 
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pattern for 
suckle pattern 


subsequent buckle 


Ripples and 


Ri g ht 


short evlinders, transition cylinders, long cylinders, 
ind very long cylinders. Short cylinders tend to be 
have as flat plate columns (cylinders with an infinite 
radius), and buckling occurs in a sinusoidal wave charac 
teristic of flat plate behavior. In contrast, the surface 
of a long evlinder buckles into a characteristic diamond 
pattern. In the long cylinder range, the critical buck 
ling stress does not depend on the length of the cylinder 
or on the end conditions. The transition cylinder be 
havior occurs between the short and long cylinder 
the effects 


length and of end conditions become important, and 


ranges. In this intermediate range, of 
buckling probably occurs in some mixed displacement 
pattern. The very long cylinder buckles by overall 
column instability without local buckling of the cylin 
der wall. The behavior of very long cylinders is ade- 
quately described elsewhere and will not be considered 
further here. 

Donnell’s eighth-order differential equation’ can be 
used to provide a small-deflection theoretical solution 


LO 


the behavior of unpressurized cylinders from the 
short cylinder to the long cylinder range. A solution 
of Donnell’'s equation by Batdorf'* gives the critical 
stress in terms of the buckling coefficient A,, which 


lor simply supported edges is defined by the equation 
\ m> + 37)? mm?) + [122Z2m? 24(m? + 8B?) 


where 








load 


a pressurized cylinder Left) Ripples at 1,500-Ib 


maximum load = 1,590 Ibs 


number of half waves in the longitudinal direc 


m 
tion 

3 L/x 

r half wavelength of buckles in circumferential 


direction 


The compressive buckling stress is computed from 


where 7 is equal to | for elastic buckling. Minimiza 


tion of Eq. (1) with respect to the parameter 


m- + B- 


TABLE 2 


Summary of Other Experimental Investigations 


Geometrical Parameters 


No.of | 
Reference Material Specimens - L 
r t r 5 

Balleretedt and Wagner Brase is 3.94, 7.87 455 4166 0.26 2.00 o 
Bridget Brass and Stee! 16 0.94, 1.68 238 919 0.70 2.04 
Clark and Holt 17S-T4 Alumioum hi) ° o29 0 5.08 
Donnell Brase and Steel 40 0.94 - 2.84 ss 469 31.83 
Fung and Sechier 24S-H Aluminum eo 5 $63 750 6.29 re} 
Kanemitsu and Nojima Stee! 86 6.38 700 304 o« 42 © 
Le, Crate, Schwarts 245-T Aluminum ? 15.00 602 6 c 5 
Lundquist Duralumis “as 5 5.00 340 4 3.00 
Robertson Stee! ° 2.50 32 6 $10 0.39 s 0 
Wilson and Newmark Stee ll 1.94 - 6.64 c 767 35 02 e 
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gives the critical buckling coefficient for long cylinders 


in the following form: 
n°)Z = 


a. = (wo 0.7022 3) 


where 2: = (rv 1 — ee 


Substitution of Eq. (3) for the buckling coefficient into 


Eq. (2) reduces to the classical equation 


0 CE(t/r) }) 
where 


C = 1/V 3(1 — p?) 0.605 for pu 0.5 
The buckling coefficient for simply supported end con 


ditions in the transition cylinder range can be deter- 
mined by substituting the values of 6B = O 


A similar solution for cylinders 


limiting 
and m = 1 in Eq. (1). 
with fixed edges deviates from the solution for simply 


T 
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supported cylinders only in the flat plate and transitj 
ranges. 

The theoretical value of the buckling coefficient k 
is given as a function of the geometrical parameter 7 
in Figs. 5-8. Theoretically, the short cylinder rang 
would occur at Z = 0. The values of A,, 1.12 an 
1.12, at Z = 1 correspond very closely to the buckling 


coefficients of simply supported and fixed-ended plat; 
columns. 
the 45 


curve 


The long cylinder behavior is represent 
by sloping straight-line portion of the curvy 
The 
ranges is the transition cylinder range 


connecting the short and the long ecvlinde 


General Analysis Procedure 


The analysis of unpressurized cylinders has consist 


of the following steps: (1) establishment of the per 


tinent parameters, (2) plotting of test data in terms oj 
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these parameters, (3) establishment of the relationship 
between the critical buckling coeflicient and the rf 


and (4) performance of a statistical analysis on 


ratio 
the test data to determine design curves. The reason- 
ing behind the choice of parameters has already been 
given 

Shown in Figs. 5-S are 252 experimental values of the 
critical buckling coefficient from the tests of twelve 
The discrepancy between theory and 
the 


investigations. 
test is clearly shown in these Figures. However 
general trend of the test points follows the theoretical 
shape of the curve quite well. Batdort, Schildcrout, 
and Stein'? have proposed semiempirical curves to ac- 
count for this discrepancy between test and theory. 
Curves at the theoretical slope were faired through the 
test data in the long cylinder range for several values of 

fon the assumption that the 7 f ratio is indicative of 
the initial imperfections The analysis employed in the 
present paper is a refinement of this procedure which 
includes a statistical basis for the correlation between 
the buckling coefficient and the y ¢ ratio. This is ae 
complished by replotting the buckling coefficient as a 
function of the 7 f ratio. Once the relationship be- 
tween the buckling coefficient and r ¢ is determined, 
the statistical analysis of step (4) is performed, and the 
design curves are drawn in terms of the theoretical 
parameters, A, and Z. 


Buckling Coefficient K. as a Function of Z 


Examination of the experimental data in Figs. 5-S 
indicates that the test points from the various investi- 
gations generally overlap, although there is consider- 
The test data have been plotted in four 
the 7 these Figures. 


Comparison of the data indicates a general increase in 


able scatter. 
increasing ranges of ¢ ratio in 
the deviation of the test values from the theoretical 


curve with an increase in the 7 ¢ range. Because of 
the large scatter and of the variation in the buckling 
coeflicient with Z, it is difficult on these plots to estab 
lish a relationship between A, and rf or to evaluate 
the significance of the scatter. It has, therefore, been 
necessary to replot the test data as a function of r 

To do this, the long evlinder test data have been 


value of the geometrical 


only. 
cross-plotted at a constant 
parameter Z. 

In the long cylinder range, the theoretical curve is a 
straight line at a 45° slope on a plot of log A, as a fune 
tion of log Z. 


cylinder range the variation of A 


It has been assumed that in the long 
with Z for each r ¢ 


ratio is given by a line parallel to, but displaced from, 


the theoretical straight line. The amount by which 


i line of constant r f is displaced from the theoretical 
curve reflects the effects of initial imperfections which 
This 


assumption will be shown to result in a family of parallel 


are assumed to be a function of the 7 ¢ ratio. 
lines for inereasing r, f. 

Each test point in Figs. 5-S can be used to establish 
a Straight line at the theoretical slope, and the test 


value of A. can be projected along this line to some 
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Compressive buckling stress coefficients as a fun 


arbitrarily chosen value of 7. A cross plot of the pro 


jected values of A, vs. r ¢ can then be made at this 
arbitrary value of Z. In reality, each line through a 
test point at a slope of 45° represents a value of the 
coefficient C in Eq. (4) different from the theoretical 


Therefore, the cross plot could also be 
with 


value of 0.605 
used to evaluate the variation of the coefficient C 


the r ¢ ratio. 


Relationship Between the Buckling Coefficient and r/t 


In Fig. 9, the projections of the experimental buckling 
coefficients are plotted as a function of the 7 ¢ ratio for 
10,000. The test points 
10,000 along 


the arbitrary value of Z 
have been projected to a value of Z 
lines parallel to the theoretical curve in the long cylin 
der range. Two ordinate scales are shown on the graph 
The A, scale is the projected value of the buckling coef 
ficient from The 
gives the equivalent values of the buckling coefficient C 


Donnell’s equation. second scale 
in the classical equation. 


Also 
classic value of 0.605 for the buckling coefficient C 


shown in Fig. 9 is a curve representing the 
5 


The large discrepancy between small-deflection theory 


and test is again clearly apparent. Comparison of the 


test results indicates that the maximum scatter 1s about 
390 per cent. 


In order to establish the variation of the buckling 


coefficient with the 7 ¢ ratio, the test results have been 


replotted in Fig. 10. In this Figure, the log of the pro 


jected buckling coefficient 4, has been plotted as a 


function of the r f¢ ratio. In addition to more clearly 


defining the relationship between A, and ? the semi 
log plot has the advantage of producing a more uniform 
scatter the The 


statistical analysis is more readily performed if the 


band for range of r ¢ considered. 


distribution of test points is uniform 

A best-fit curve is shown on the log 4. vs. r ¢ graph. 
The shape of this curve was selected by the investiga 
tors as being representative of the trend of the experi 


mental data. The vertical location of the curve with 








594 JOURNAL OF THE 


AERONAUTICAL 





BALLERSTEDT & WAGNER * DONNELL = LO, CRATE & SCHWART/ 
* BRIDGET 2 FUNG & SECHLER * LUNDQUIST 
+ BUGHY (NAA = NAA « ROBERTSON 
* CLARK @ HOLT @ KANEMITSU & NOJIMA « WILSON & NEWMARK 
4 
e 
. ~ THEORETICAL CURVE ‘a 
ww C0605) BEST FIT CURVE — 
“Sark 
peek = t : 
2 we a nia Oo ° "le Ge 
CC ~g._2 + a 
u ~~a, *g° eo = - el " , 
4 t 4 ya! eo? t ¥ = 
a 3 a) —_ : — 
’ e <—s SS 
q as a . 
aR i ceo oe * 
M 
. 90% PROBABILITY CURVE 
99% PROBABILITY CURVE 
50 1000 500 2000 2500 3000 3500 4006 450¢ 500 
wt 
Fic. 10. Compressive buckling stress coefficients as a function 


of y ¢ (semilogarithmic grid 


respect to the test data was achieved by satisfying the 


equation 


aC V 0) s 
i 
where 
V log A. at Z 10,000 for test points 
v log A, at Z 10,000 for points from the best 
fit curve at r fof test point 7 
If Eq. (5) is satisfied, it can be shown that the follow 


ing relationship is also satisfied: 
a minimum (G) 


Thus, by satisfying Eq. (5), the assumed curve was in 


the least square position. 


Statistical Analysis and Development of Design Curves 


After a best-fit curve had been established with the 
sum of the squares of the deviations a minimum, a 
IBM 701 


digital computer was used to solve the equations of the 


statistical analysis was performed. An 


statistical analysis in the following form: 


Vo = Va Co (7) 
where 
Vs log A, at Z = 10,000 for a specified probability 
J V (vi — ¥.)? (N — 1) = the standard deviation 
NV. = number of test points 
on probability coefficient 


The solution of Eq. (7) for the 222 test points gave a 
value of 0.138 for the standard deviation of the log of 
K,.. The value of C, = 1.396 was selected from a small 
sample statistical table*’ for a 90 per cent probability 
The 90 per cent probability design 


design curve. 


curve is drawn in Fig. 8. This curve agrees well with, 
and generally parallels, the lower bound of the test 
data. 


on this graph. 


The 99 per cent probability curve is also shown 
The relatively large difference between 
the 90 per cent and the 99 per cent probability curves 
is typical of a statistical analysis of test data with a 


large scatter. 
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Final 90 per cent probability design curves were ob 


tained in the following manner: For a_ particular 


10,000 was 


read from the 90 per cent probability curve in Fig. 10 


value of r ¢t, the design value of A, at Z = 


This value was plotted on a graph of A, vs. Z, and 
straight line at the theoretical slope of 45° was drawn 
through the point. At low values of the Z parameter 
the buckling coeflicient approaches the wide plate 
column value which depends on the end conditions 
For design purposes, values of | and 4 have been as 
sumed for simply supported and fixed ends, respe 
tively, at a value of Z = 1. A transition curve has 
been fitted between the long cylinder range defined by 
the straight line and the flat plate value at 7 

The procedure suggested in reference 21 has been 
used to construct the design curve in the transition 
range. In the authors’ investigation of the compres 
sive buckling of curved panels, sufficient data were 
available in the transition region to indicate that the 
method would give a reasonable design curve. Al] 
though there were not enough data to check the transi 
tion curve for cy linders, the procedure was adopted on 
the basis of its success in the case of curved panels 


Basically, this transition curve 1s defined by 


A Ky (12Z?/#*K 51) \ 


where A, the simply supported or fixed-end wide 
plate column buckling coefficient. This equation can 
be made to coincide with the flat plate values at 7 

and is tangent to the design curve in the long cylinder 
range. The complete design curves for the range of 
r t likely to be encountered in aircraft structures are 
shown in Fig. 11. A comparison with the theoretical 
curve indicates the large difference between theory and 
the 90 per cent probability design curves. For in 
stance, the design curve for an r ¢ of 3,000 is only about 
12! . per cent of the classical theoretical values in the 
long column range. The 90 per cent probability curves 
can be compared with the experimental data in Figs 


5-8. In general, these statistically defined design 
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curves are significantly lower than previously available 


lesign curves. 


Statistical Design Criteria 


The load at which buckling of a thin monocoque 
shell occurs can be considered as ultimate since no 
further load-carrying capacity remains. Because the 


final design curve of this study is defined in terms of a 
statistical probability of failure, it is incumbent on the 
user of these data to decide what probability is required 
by the design criteria demands of his particular struc 


ture In some ways, this is a new concept for the 


inalyst, since in the past many accepted analytical 


procedures could not be conveniently related to failure 


probability. 

rhe structures engineer usually has two design loads 
to consider limit load and ultimate load. In general, 
limit load is defined as the maximum expected operat 
ing load. Ultimate load is usually defined as the limit 
load times a factor of safety which accounts for uncer 
Ultimate load 


may normally be experienced by an aircraft only in a 


tainties of load, analysis, material, ete. 


static test. Traditionally, piloted aircraft have re 
quired larger factors of safety than unpiloted missiles. 

lo illustrate the placement of design curves and to 
the 


philosophy, 90 and 99 per cent probability curves have 


discuss weight-saving possibilities of a failure 


been plotted in Fig. 10. These probabilities are iden 
tical to those used for determining the A and B ma 
terial properties in ANC-5. 

\s shown in Fig. 10, the 99 per cent probability curve 
is about 30 per cent lower than the 90 per cent probabil 
ity curve. With reference to Eq. (7) it can be seen 
that this would be characteristic of data with consider 
able scatter which causes a large standard deviation. 

If the 99 per cent probability curve is used for ulti 
mate load, there is a slim 1 per cent probability of fail 
ure at ultimate and practically no possibility of failure 
it limit load. If the 90 per cent probability curve 1s 
selected, there is a 10 per cent probability of failure at 
ultimate load and a probability of failure at limit load 
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which would probably be less than 2 per cent depending 
| ; | 


on the ultimate to limit ratio \ccompanying the 10 


per cent risk, a significant weight saving in the shell is 


achieved. Since present missile reliability philoso 


phies are being established on a statistical basis, this 


type of analysis offers attractive weight-saving possi 


bilities. 


ANALYSIS OF PRESSURIZED CYLINDER TEST [A1 
Theoretical Considerations 
Experimental investigations have shown that the 


critical buckling stress is significantly increased by 


internal pressure. Application of the small-deflection 
theory to the compressive stability of cylinders with 
been unsuccessful. However 


internal has 


Lo, Crate, and Schwartz® have analyzed the problem 


pressure 


of long pressurized cylinders using an extension of the 


large-deflection theory of von Karman and Tsien 


Plotting their results in terms of the nondimensional 


parameters (P E)(r t)* and (¢., /)(r t), they tound 
that the buckling coefficient C increased from the 
Tsien® value of 0.375 at zero pressure to the maximum 
classical value of 0.605 at (p E)(r 1) 0.169. The 


large-deflection theoretical curve is shown in Fig. 12 
Also shown are the experimental values of Lo, Crate 
and Schwartz with the experimental data obtained in 
the present investigation and that of reference I5. It 
is apparent that, as for unpressurized cylinders, large 
discrepancies exist between the theoretically predicted 
and the experimentally determined buckling stresses 


For example, in the region of maximum critical stress 


the experimental values of C vary from about 0.5 to 
0.585, whereas theory predicts 0.605. Lo, Crate, and 
Schwartz suggested that better correlation with ex 


perimental results could be obtained if the increment in 
the 
plotted as a 


stress parameter Ao.,'f t were 
the 


13. The increase in 


buckling 


function of pressure parameter as 


shown in Fig. the critical stress 
Ag; 


beneficial effect of internal pressure. 


due to internal pressure directly represents the 
The total critical 


stress is obtained by adding the critical stress for the 
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unpressurized cylinder to the increase in the critical 
stress due to pressure. 

The theoretical curve and the test results are shown 
in Fig. 13 in terms of the incremental stress parameter. 
In order to plot these data it was first necessary to de- 
termine the unpressurized critical buckling stress for 
The 90 per cent probability design 
As can 


each 
curves of Fig. |! were used for this purpose. 
13, the general trend of the experi 


specimen. 


be observed in Fig. 
mental data agrees quite well with the theoretical 
curve. It was this agreement between test and theory 
which suggested the feasibility of a semiempirical analy 


sis of the data using these parameters. 


Development of Design Curve 


The analysis of the pressurized cylinder data was 
achieved by first selecting a best-fit curve which is 
shown in Fig. 15 together with the experimental data. 
This curve has been drawn to approach an asymptote 
at large values of the pressure parameter and to de 
crease at low values with essentially the same slope as 
the theoretical The between the 
best-fit curve and the experimental data is apparent 
The best-fit curve has been drawn so 


curve. agreement 
from Fig. 13. 
that the sum of the squares of the deviations of the ex 
perimental data from the curve is a minimum. 

A statistical analysis of the test data was performed 
to establish the 90 per cent probability design curve of 
Chis design curve is to be used in conjunction 
The 


curve approximates a lower bound to the test data. 


Fig. 13. 


with the unpressurized design curves of Fig. 11. 


The standard deviation of the log of the stress parameter 

O.116. 
the proposed design curve at large values of the pressure 
For small values of the pressure parameter, 


The theoretical curve nearly coincides with 
parameter 


the design curve is appreciably lower than the theoreti 


cal curve. 


CONCLUSIONS 


The axial compressive buckling of thin-walled circu 
The 


oretical parameters and shapes of curves were used and 


lar cvlinders has been analyzed semiempirically. 
the test data were fitted statistically. The basic ap 
proach to this method is that the initial imperfections 
of the test evlinders were equivalent to those which 
would be encountered in aircraft structures. Design 
curves based on a 90 per cent probability criterion have 
been presented. For unpressurized cylinders, the 
critical buckling stress coefficient for short and long 
cylinders is given as a function of the geometrical param- 
eters Z and r ¢. For long, pressurized cylinders, the 
increase in the critical stress due to pressure is given as 
a function of pressure and geometry. In general, these 
statistically defined curves are significantly lower than 


TIC 
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previously available design curves. It is proposed that 
the use of statistically defined buckling curves can fy 
correlated with missile reliability philosophies 
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Effects of Leading-Edge Separation on Thin 
Wings in Two-Dimensional 
Incompressible Flow’ 


L. F. CRABTREE* 
Royal Aircraft stablishment, Farnborough, Eingland 


SUMMARY 


| flow over an airfoil with a laminar separation bubble near 

e leading edge is treated. There are thus two attachment lines 

ind two separation lines, in contrast with the viscous flow over 
irfoil with separation only near the trailing edge and with 

l ir the leading edge 

A physical picture of such bubble flows is constructed, and a 


mis sought by splitting up the problem into three parts 





concern the viscous flow inside the bubble, the turbulent 
mixing region of reattachment, and the inviscid external flow 
Little is as vet known about the internal flow of a bubble, but 
the other two problems are treated in some detail 


paid to the problem of transition in the 


Schubauer and Klebanoff’s condition for 


Attention is also 
separated shear layer 
srowth of turbulent spots in a boundary layer is employed 


the g 

to explain the suecess achieved by Owen's criterion. This was 
suggested as a means of distinguishing between the two types of 
bubbles (‘long” and ‘‘short’’) commonly observed on thin air 


foils at incidence, 
rhe picture of laminar separation bubbles thus obtained offers 
1 explanation of the differences in behavior of long and short 
, and a hypothesis concerning the leading-edge stall or 
brupt change from a short to a long bubble is put forward. 
the problem of predicting maximum lift is discussed 
while this is still not possible in any specific case, the scale effect 
i the stalling characteristics of thin and moderately thin airfoils 
is explained in some detail 
Phe whole paper is in the nature of a summary report, how 
Royal Air 


craft Establishment on the problem of regions of separated flow 


ver incomplete, of recent work done mainly at the 


ensional airfoils at low speed 


INTRODUCTION—-ViIscousS FLOW WITH SEPARATION 


ONLY AT THE TRAILING EDGE 


bien rHEORY of viscous flow over an airfoil with 
separation only at the trailing edge is now fairly 
well established. It will be recapitulated here in order 
to emphasize the comparison with the flow over an air 
foil including regions of separated flow on the surface. 
In the former flow there is one attachment line near the 
leading edge and one separation line near the trailing 
sketched The 
viscous region around the airfoil in the form of a 


edge, as in Fig. 1(a) presence of a 
boundary layer and wake does not change this type of 
flow but affects primarily the shape of the displace 
ment surface and the value of the circulation. 
Presented at the Aerodynamics—-I Session, Twenty-Fifth 
Annual Meeting, IAS, New York, January 28-31, 1957 
rhe author is deeply grateful to his colleagues at the Royal 
Aircraft Establishment, and in particular Dr. D. Kiichemann, 
lor their advice and encouragement in the preparation of this 


Senior Scientific Officer, Ministry of Supply 


The circulation about an airfoil may be predicted 


from detailed calculations of the boundary layer and 


wake; this was first done successfully by Preston 

Some improvements of this method have been effected 
recently by Spence rhe theory rests largely on 
G. I. Taylor's theorem regarding the overall circulation, 
and the lift force on the airfoil may be calculated once 
predicted and the 
Maskell® has shown 


term; 


the overall circulation has been 
circulation in the wake subtracted. 
that the 
it is given by 

The 
airfoil can be found, and from it the pitching moment 


circulation is usually a small 
Voe-Cp-a 


distribution 


wake 
for a symmetrical airfoil. 
pressure over the surface of the 
and the position of the aerodynamic center, by allow 
ing for the displacement effect of the boundary layer 
This may be done by simply adding the displacement 
thickness of the boundary layer to the actual profile 
shape--1.e., the shape of the boundary of the external 
inviscid stream is taken from the solution of the vis 
cous flow but the boundary condition is still the same, 
namely, that there is no normal velocity. An ex 
ample where all the various steps have been performed 
in detail from experimental data has been given by 
Brebner and Bagley. 

In principle the drag force due to viscosity can also 
be determined. A convenient practical method is due 
to Squire and Young,’ and simpler developments of this 
theory have been proposed by Hoit® and Tetervin.’ 
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Iwo types of viscous flow over an airfoil 
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Fic. 2. Laminar separation bubble on an airfoil 


Frow Witrn A BuBBLE DUE TO LEADING-EDGE 
SEPARATION 


The type of flow considered so far is no longer ade 


quate when separation occurs elsewhere than only near 


the trailing edge —as, for instance, near the leading edge 
on some airfoils near maximum lift or on thin airfoils 
even at moderate or low values of lift. This is the 
case we wish to consider here. The flow with only 
trailing-edge separation is then no longer possible; 
it would imply a high suction peak near the leading 
edge followed by a severe adverse pressure gradient 
which would lead to a separation of the laminar 
boundary layer there. The viscous region now loses 
the character of a boundary layer since, as Maskell 
points out,* the so-called singular separation along a 
line normal to the stream must lead to the formation of 
a bubble in this case. The bubble includes air which 
is separated from the mainstream and is carried along 
with the body. We shall consider the case where the 
bubble length is smaller than the wing chord—4.e., 
the flow reattaches to the surface. Thus, as sketched 
in Fig. l(b), there are two attachment lines, A, and 
two separation lines, S, in direct contrast to Fig. I(a). 

We now consider how the flow in and around such a 
bubble is likely to look, and Fig. 2 has been constructed 
as follows. We assume a closed bubble so that a ‘‘zero”’ 
streamline leaves the surface at S and joins it again at 
A. Further streamlines from infinity upstream stay 
outside this boundary. The boundary layer is of 
finite thickness at S, and since there is no longer a 
boundary condition requiring zero tangential velocity 
along the zero streamline (as there is on the solid sur 
face) the ‘“‘shear layer’? must diffuse and extend inside 
the zero streamline. The air is thus set into motion 
inside the bubble and for a closed bubble continuity 
demands a circulatory flow inside the bubble. 

The flow over the front part of the bubble is entirely 
consistent with laminar flow in the viscous region. 
If the bubble reattaches, however, there must be a 
considerable widening of the streamlines over the rear 
end. This cannot be imagined as a laminar process at 
the Reynolds Numbers we are concerned with, but is 
more akin to turbulent mixing as it occurs in, say, a 
pipe with a sudden enlargement of cross section. We 
are, therefore, led to assume that there is transition 
to a turbulent state and that this is likely to occur near 
station (1) or slightly further upstream. We also 
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arrive at velocity profiles as shown in Fig. 2; since ther 
is zero skin friction at station (2) the profile there wil] 
show the vertical tangent characteristic of a separation 
profile. The velocity profile at station (2), however, is 
not necessarily identical with the usual separation pr 

file of a fully developed turbulent boundary layer 

With this picture of the flow we note that there ar 
likely to be large pressure changes (in fact, a pressur 
recovery) between P; and A, but a nearly constant pres 
sure between S and P). 

We note also, in searching for a solution, that this js 
not like an ordinary boundary-layer problem. It may 
be possible, however, to retain some of the assumptions 
of boundary-layer theory; for example, that inside the 
viscous region there is no variation of static pressure 
with distance normal to the surface. It may also be 
possible to proceed on the assumption that the problem 
can be subdivided into three partial problems namely, 
solutions for the three main regions 1) inside the 
bubble, within the zero streamline; (2) the turbulent 
mixing region, the major part of which may possibly 
be limited to P, P,P; A; and (3) the external inviscid 
stream along the displacement surface, the pressur 
distribution along which may be taken as approxi 
mately that along the wall. 

Little work has been done as yet on the internal 
flow of the bubble. 
and Batchelor’ on = similar problems assuming a 


However, the work of Squire 


“core” of constant vorticity surrounded by a layer in 
which boundary-layer assumptions can be made may 
be applicable and fruitful. This work should be ex 
tended to include a consideration of the (assumed 
laminar surrounding layer so as to investigate its sta 
bility and determine the transition point. 

The mixing region will be considered later; this 
may be used to estimate the order of the pressure re 
covery and hence the level of suction in the bubble, 
and also the energy transport across the zero stream 
line P; A needed to maintain the flow in the bubble. 
Strictly speaking, of course, such quantities can only 
be found from a complete solution of the entire flow 
field. The separate solutions for the three regions can 
not be unique in themselves, and it is onlv by match 
ing all three that a unique solution can be found. A 
rough answer can be obtained from any one region 
though, provided its interaction with the other two 1s 
adequately represented. 

The problem of the external stream will also be 
considered more fully later. This again differs from 
the earlier problem of airfoils with ordinary boundary 
layers because the shape of the displacement surtact 
is not known to begin with. Following Maskell, 
we propose a different course: use is made of two 
experimental observations (a) that the pressure 1s 
constant over the forward part of the bubble, and (b 
that the displacement thickness of the viscous region 1s 
of the order of the bubble thickness downstream of the 
bubble. These lead to the assumptions that the region 
(a) extends back to the maximum height of the bubble, 
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that the variation of the displacement thickness aft 


ind that 
{this point is negligible. 
Finally, as emphasized above, the three solutions 


% joined together and matched. For instance 


must | 
the basic solution for the external stream is not unique 
in that the level of suction in the bubble is undeter 


Chis will be obtained by matching the pres 


mined 

sure recovery With that in the turbulent mixing region. 
Cure EFFECT OF TRANSITION IN THE SHEAR LAYER 
So far no characteristic dimension of the flow has 

heen fixed. The size of the bubble—1.e., its length in 


relation to the airfoil chord~—must depend very much 
m whether transition occurs a long way downstream 
{ separation, S, or almost instantaneously. In tests 
fa flat plate held normal to the stream Flachsbart 


found that transition in the free shear layers behind the 


plate jumped forward suddenly at some critical Rey 


nolds Number. This is also related to the critical 
Reynolds Number for transition in the wake behind a 
evlinder observed by Linke. We may, therefore, 
expect a discontinuous change in the length of path 
from S to P; (sav) needed to produce transition to a 


turbulent state, and, therefore, either a very short 
bubble or a fairly long one. It may also be expected 
that the sudden switch-over depends on a Reynolds 
Number related to conditions at separation, for ex 
ample that based on boundary-layer displacement 
thickness at S, (1-61/v)s (R51) s 

rhis is confirmed by experiment. Owen and Klan 
fer!‘ have analyzed tests in which bubbles were observed 
with this in mind, and found a relation betwcen bubble 
length and (R5 It appcars that there are two dis 
tinct tvp-s of bubble. 

A short bubble for which (R35) 5 1s always greater 
than about 500. Its length is of the order of 1 per 
cent of the chord or less [say 100 times the boundary 
layer displacement thickness at separation (6;)s| and 
the bubble contracts as the angle of incidence is in- 
creased 

2) A long bubble for which (R6)s is always less 
than about 500, and whose length is of the order of 2 
or 3 per cent of the chord on formation at low incidences 
the ratio of bubble length to (6;)s5 being of the order of 
10? or more The long bubble grows rapidly with 
increasing angle of incidence until the separated layer 
fails to reattach to the airfoil surface and the bubble 
extends out into the wake. 

Owen's criterion has been confirmed by the analy- 
sis’:'® of more recent tests, but it appears that a short 
bubble does not always break down into a long bubble 
at the same value of (Rs,)s. Since the breakdown also 
gives rise to a marked change in the external flow field 
this is not altogether surprising. In a sense the Owen 
criterion by distinguishing between the two types 
ol bubble —does tell one that the breakdown will be 
abrupt. We must have one kind or the other, with no 
intermediate state. The success achieved in distin- 


guishing between the two types of bubble can be more 
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completely explained, however, on the basis of Schu 
bauer and Klebanoff'’s condition” for the growth of 
turbulent spots. In their experiments on the boundary 
layer on a flat plate as part of an investigation into the 
mechanism of transition they found that turbulent 
spots grew only very slowly if Rs, was less than about 
150, and a much faster rate of growth was obtained 
when 5, exceeded this value. Now MeGregor" has 
found that the displacement thickness of the laminar 
shear layer of a bubble grows only very slowly from 
separation up to the onset of transition. It is clear, 
therefore, that, if (R;)s is less than about 450-500, a 
long path from 5 to P; (say) is required for transition, 
whereas if (Rs) 5 > 450-500 transition may occur quite 
rapidly. In the first case a long bubble would be 
formed, and in the second a short bubble. 

The actual mechanism of the abrupt change from a 
short to a long bubble (of which the leading-edge stall 


is a particular example) will be discussed later 


THE TURBULENT MIXING REGION 


A simplified model of a type of flow with a bubble on 
1 two-dimensional airfoil has been considered by Nor 
bury and Crabtree,'’ and particular attention was 
given to the reattachment process. They assumed that 
the principal mixing process and the corresponding 
pressure recovery occur in the region P)P2P;A of Fig. 2 
be 


linear, while that at A was assumed to be of the char 


The velocity profile across P,P. was assumed t 


acteristic separation type with a vertical tangent at the 
wall, the profile being represented by a simple power 
law. It was also assumed that fluid passing through 
P;P. at station (1) is diffused to form a boundary 
layer of depth AP; at station (2) while the velocity ot 
the flow at the upper boundary of the mixing region 
falls from 1) to I without loss of energy. To simplify 
the analysis the following further assumptions were 
made. 

(1) P2P; is a streamline. Further, no fluid enters or 
leaves the region across BP; on the average, from the 
definition of the bubble boundary. 

(2) The curvature of the airfoil surface SBA may be 
neglected. 

3) P2P; is parallel to BA. 

Applying the momentum theorem and the equation 
of continuity, an expression was obtained for the co 
efficient of pressure rise, ¢, over the rear part of the 


bubble, where 
o = (P. — P,)/(1/2)a V;? = (Cp Cx | c. 


It was found that the value of ¢ varied only slightly with 
the ratio of thicknesses 6,/6, (where 6, P,P» and 
64 = AP 
For instance, if the latter 1s given by 


and with the velocity profile across P,P 


vy Vy a+ l a b On 


and the profile at A is given by 
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in the limiting case » = 1, o varies from 0.37 to 0.45 
as the constant a varies from 1.0 (a uniform profile) to 
0.5. 

A further way of predicting a reasonable value of o 
is suggested by a consideration of the rate of spreading 
of the turbulent mixing region between two parallel 
streams of different velocities. This problem has been 
analyzed by Kuethe,*’ and from his results it may be sur- 
mised that when the pressure rise through turbulent 
mixing ceases (roughly with the return of the zero 
streamline to the wall) the depth of the mixing region, 
AP, in Fig. 2, will be about three times as great as the 
original depth of flow involved in the mixing (P,P. in 
Fig. 2). This yields a value of o of 0.52, which is of 
the same order as before. Also the fact that the tur 
bulence spreads at a much greater angle into the more 
slowly moving stream supports the assumption that 
PP; is parallel to BA. 

The analysis of pressure distributions obtained by 
Dimmock*! in tests of a jet-flapped airfoil where the 
long bubble is fairly well defined shows that the value 
of ¢@ is in fact roughly constant at 0.55. 

In view of the assumptions made in the analysis this 
value of o may be considered to be an upper limit. 
The pressure recovery may be less than this in some 
cases, but it cannot be much greater. It may be re- 
called that the pressure recovery due to complete mix 
ing after a sudden enlargement in a pipe yields a maxi- 
mum possible value for o of 0.5 (for an area ratio of 
0.5). 

In the particular case of a short bubble, Crabtree!® 
has shown by analysis of the experimental results in 
references 1S and 22 that the value of o gradually rises 
with increase of incidence (or decrease in Reynolds 
Number at constant incidence). When a limiting value 
of about 0.35 is reached the short bubble suddenly 
“bursts” and the flow pattern breaks down into that 
associated with a long bubble. This has been pro 
posed as part of the mechanism of the leading-edge 


stall, and will be more fully discussed later. 


THE EXTERNAL INVISCID FLow 


The problem of the external flow has been treated 
successfully by Maskell'! for the simple case of a two 
dinensional flat plate. The presence of a long bubble 
necessarily implies a thick boundary layer behind the 
bubble and a thick wake, neither of which can be neg 
lected in a treatment of the external flow problem. 
Maskell’s displacement boundary, therefore, consists of 
a constant pressure region-——assumed to extend up to 
the maximum thickness of the bubble— followed by a 
constant displacement thickness up to the trailing edge, 
and a wake of nearly constant thickness which leaves 
the trailing edge smoothly and curves back to the main 
stream direction at infinity downstream. 

The model gives rise to a rapid pressure rise just 
downstream of the constant pressure region, after 
which the pressure distribution steadily approaches 
that in the absence of a bubble i.e., for the flow with 
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one separation line and one attachment line s‘sctched 
in Fig. l(a). 

Clearly a long bubble could not exist with the sany 
value of peak suction (C, as in the flow of Fig 
at the same incidence: the required pressure ris 
could not be supplied by turbulent mixing. Thx 
suction, therefore, collapses. Now let us consider th 


} 
f 


effect of an increase in incidence for this type of fle 
If the level of suction in the bubble remains unchanged 
then in the above model the bubble length increas:s 
as sketched in Fig. 3. This also fits with the hypothesis 
of a constant pressure recovery o in a long bubbl 
In fact, the level of suction usually falls with increasing 
incidence, which means that the bubble expands « 
more rapidly. 

The effect of incidence on a short bubble is represented 
schematically in Fig. 4. In this case boundary-laver 
calculations indicate that separation moves forward 
toward the leading edge. This can be seen, for in 
stance, by application of the simple criterion that the 
ratio of velocity at separation to the maximum velocity, 
I’s/Vn, 1s about 0.94 in a typical case of this sort. 
The levcl of suction in the bubble is, therefore, higher 
at a higher incidence. This implies a more highly 
curved flow to maintain the increased suction, and, 
thercfore, the bubble is shorter. Here, turbulcnt 
mixing can supply enough pressure recovery to bring 
the pressure back to that which would exist in the 
absence of the bubble. Therefore, while the bubble 
shortens as incidence incr:ases the pressure distribu 
tion is much the same as in attached flow, except for 
a very small region. The bubble is only a means of 
producing a turbulent boundary layer and making 
possible a flow which is, practically, ‘‘attached” near 
the leading edge. 

These arguments are largely confirmed by experi 
ment. That the short bubble has only a very small 
effect on the pressure distribution is shown by the 
comparison in Fig. 5 where the experi:nental pressure 
distribution is taken from tests on a NACA 65-009 sec 
tion by McCullough and Gault.** Consequently the 
values of lift, drag, and pitching moment are likewis 
little affected by the presence of a short bubble 

A preliminary comparison between Maskell’s thcory 
of a long-bubble flow and experimental results for a 
NACA 64A-006 section in Fig. 6 shows that the pres 
sure distribution is well represented provided the thick 
ness of the airfoil is taken into account. This was donc 
by the simple expedient of subtracting the measured 
pressure distribution at zero lift. As stated above the 
value of the suction in the bubble, on which the shape 
depends, is left open in the theory. A unique solution 
can be obtained by applying the condition obtained 
from an analysis of the turbulent mixing region, that 


the pressure recovery coefficient o is equal to 0.5¢ 
This value of « may be used to choose the physically 
possible solution from Maskell’s infinite number of 
mathematical solutions such that the pressure <ls- 
tribution downstream of the bubble appreaches that 
calculated for attached flow at the same lift coeffi- 
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cient. We have assumed that in the case of Fig. 6 
the turbulent mixing is practically over by about 45 
per cent chord where the pressure and loading dis- 
tributions with bubble approach those for attached 
flow. Thus C,, = — 0.35 and C, (C,, — «)/(1 - 
a) = — 1.1 (for ¢ = 0.35) which leads to the solution 
drawn. The difficulty always lies in defining C,,—1.e., 
the location of reattachment. The example does 
show, however, that there is some prospect of develop- 
ing and combining these methods into one of practical 
utility. 

The long bubble, of course, has a great effect on the 
pressure distribution, and some further examples are 
given in Fig. 7, from tests by Butler.*4 The large 
drag increase associated with a long bubble obviously 
stems from the breakdown of the high suction peak 
near the leading edge, and the shift of the center of 
pressure is linked with the gradual lengthening of the 
bubble as the incidence is increased. Note that the 
trailing-edge pressure is appreciably reduced when 
the length of the constant pressure part of the bubble 
exceeds about half the chord. 

This brings us to one of the limitations of our pic 
ture of bubble flows: if reattachment (A in Fig. 2) 
approaches the trailing edge, then the static pressure 
there falls since not a// of the mixing occurs upstream o* 
the trailing edge. Consequently, (1) there is more 
circulation in the wake, and, therefore, less lift, and 
(2) the circulation is not “‘fixed”’ so that fluctuations 
are much more likely to occur. This leads to buf- 
feting a problem of great current interest, but which 
is unfortunately outside the scope of this paper. 


MAXIMUM LIFT. AND STALLING CHARACTERISTICS 


The influence of the type of bubble formation on the 
aerodynamic characteristics of an airfoil section is now 
fairly clear. In particular, the effect on stalling char 
acteristics (defined here as the behavior at maximum 
lift) has been described in a classic report by McCul 
lough and Gault.** 

Altogether there are three main types of stall: (a) 
trailing-edge stall, with the separation point of the 
turbulent boundary layer moving forward from the 
trailing-edge as the incidence increases; (b) leading 
edge stall, caused by an abrupt separation of the flow 
near the nose without subsequent reattachment —.e., 
short bubble ‘“‘bursting’’; and (c) thin-airfoil stall with 
laminar separation near the leading-edge and turbulent 
reattachment at a point which moves progressively 
rearward with increasing incidence 1.e., a long bubble. 

On the basis of previous discussions a hypothesis con 
cerning the leading-edge stall has been proposed." 
Thus the sequence of events in the breakdown of a 
short bubble could be as follows. 

(1) A short bubble is no longer possible, because the 
maximum possible value of o has been reached so that 
the flow cannot reattach immediately. 

(2) The peak suction collapses, and, as a result of 
the redistribution of pressures around the nose and the 


consequent effect upon the development of the laminar 
boundary layer, the value of (Ré;)s falls to below, say, 
500. 

(3) Thus turbulent spots are unable to grow and 
the separated shear layer remains laminar for a much 
greater distance. 

(4) Eventually, however, the value of Rs, of the 
separated layer exceeds 500 so that transition ensues 
and finally turbulent mixing occurs. 

(5) Now the suction in the bubble is low enough 
and the pressure distribution is such that the bubble 
can close with a small enough value of the pressure re- 
covery factor a. The closure of the bubble, however, 
may occur behind the trailing edge. 

Thus all three regions of the bubble flow as outlined 
above have accommodated themselves again so that a 
possible flow results. The loss of lift is abrupt on the 
breakdown of a short bubble if the closure of the sub 
sequent long bubble is near or aft of the trailing edge 

We must bear in mind, however, that the condition 
for the growth of turbulent spots was formulated for 
an attached flow; the situation may be different for a 
shear layer away from the solid surface. The problem 
of the stability of laminar wakes 1s being studied by 
MecKoen” and Curle,*® and their results should later 
be applied to our particular problem. 

Long bubbles, on the other hand, must lead to an 
abrupt change in the lift force on formation. As they 
lengthen with incidence, the rate of increase of the over 
all circulation is gradually slowed down and the circu- 
lation in the wake increases, so that the lift curve de 
parts from the straight line and bends over. Long 
bubbles can have significant effects on the drag and 
pitching moment characteristics of an airfoil as well 
as on the lift, as has already been pointed out. These 
have been discussed by Kiichemann™ and Crabtree." 
There must also be differences in the actual values ol 
maximum lift coefficient C,,,,,, achieved in the different 
types of stall. Clearly, the lift force on an airfoil sec- 
tion with long bubble will always be below that on an 
airfoil with short bubble and so will their respective 
maximum lift values, as is indeed observed in experi 
ments. The prediction of the actual value of Cy,,,, 18 
very difficult, and, as yet, only general trends can be 
noted. For example Multhopp® has constructed 
curves of C,,,,, against foo5, the ratio of thickness at 
0.05 chord to the chord length, for a number of airfoil 
sections, while others have used the nose radius as an 
empirical parameter. 

We are now also in a position to discuss the possible 
scale effect on the maximum lift of airfoils experi 
encing bubble separations. Possible lift curves for a 
thin airfoil at three different Reynolds Numbers are 
sketched in Fig. 8. For the lowest value of R the 
curve (c) shows a kink at a certain (fairly low) in 
cidence, representing the formation of a long bubble 
which still reattaches upstream of the trailing edge 
An increase in R shortens the length of laminar sepa 
rated layer over the front part of the bubble so that the 
pressure-recovery coefficient o is reduced. Thus a 
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higher incidence may be attained before the maximum 
ilue of o is reached and the bubble bursts into 


pe yssibl 


the long type, which is then accompanied by a loss of 
lift. If the Reynolds Number is sufficiently high the 
change in the type of bubble formation may be post 


to such a high incidence that the resulting long 
(or even extends out 


hubble covers most of the chord 
into the wake) so that maximum lift had been at 
tained. Curve (b), therefore, illustrates a_ typical 


leading-edge stall. 
On the other hand, the thinner boundary layer ob 


| at a higher Reynolds Number means that there 


tamed < 

is less loss of lift due to viscous effects. Therefore, at a 
given incidence the suction peak is higher and the fol 
lowing adverse pressure gradient is steeper. The pres 
sure recovery o is higher on this score so that in some 
cases an increase in R may mean that the short bubble 


breaks down at a lower incidence. Although this may 


be expected to be usually a secondary effect, cases 


where it is predominant have been observed e.g., 
in \MIeCullough’s tests?* of a NACA OO06 section. 
he curve (a) representing the highest free-stream 


Number shows a typical rear-stalling be 


lee \ nolds 


havior. For such a case, it may be that FR is so high 
that transition to turbulence in the boundary layer oc 
curs ahead of laminar separation. Crabtree* has 
given as an empirical condition for this to be so that 
R35 must be greater than about 2,700. The first 


separation other than at the trailing edge is then one 


of the turbulent layer somewhere at the rear of the sec- 


tion This leads again to a bubble extending into the 


wake, and its gradual growth with angle of incidence 


causes the lift to reach a maximuin value. Also some 


thick and some cambered airfoils stall in this manner 
while still exhibiting a short bubble near the leading 
edge. The rear stall is then accelerated by the thicker 
boundary layer which results from the short leading 
edge bubble Similar ideas have been put forward in 
dependently by McGregor." 

that the curves of Fig. S 


represent the same airfoil at different Reynolds Num 


It must be remembered 


bers. It is still possible for a thin airfoil with leading 
edge stalling behavior to attain a higher C, than 
does a thicker airfoil with a rear-stalling behavior. 


In Fig. S, however, the rear-stall pattern produces the 
highest C, 

At a given incidence below the stall the difference in 
lift coeflicient C, between the long bubble flow of curve 
¢) and those at higher Reynolds Numbers may be ex 
plained, according to Maskell,” almost entirely by a 
consideration of the wake circulation. He does in fact 
find that all the main characteristics of the lift-curve in 
f a thin-airfoil stall can be explained in terms 
i.e., the break in the lift-curve and 
the subsequent loss of slope. Whatever the type of 
stall though, the flow patterns should eventually be 


the case « 


of wake circulation 


very similar for the same airfoil and so we expect all 
three lift curves of Fig. 8 to run together beyond their 
respective incidences for maximum lift. 

However, although it is impossible to calculate the 


SE 
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complete lift curve in any specific case (including the 


prediction of C,,, we can design a section on which 
laminar separation with a long bubble is avoided, up to 
a point. This can be done by specifying the basic pres 


sure distribution, and camber is a suitable means ol 


fulfilling such requirements. Boundary-layer control 
can also be effective, and either suction or blowing out 
of small holes has been suggested (see reference 16 


CONCLUDING REMARKS 


A fairly complete picture of the phenomenon of lami 
nar separation bubbles on thin airfoils at incidence is 
now available, but a number of points, such as_ the 
viscous flow inside the bubble, require further investi 
gation. The complications of the whole problem have 
been emphasized by comparison with the well-estab 
lished theory of viscous flow over an airfoil with s¢ para 
tion only near the trailing edge. 

Only the case of two-dimensional airfoils has been 
considered, so that a major problem which remains 1s 
the study of three-dimensional separations 

Of the two types of lamimar separation bubble com 
monly observed, long and short, the former 1s perhaps 
the more significant due to its large effect on the aero 
dynamic characteristics of an airfoil. Short bubbles 
must also be considered further in their capacity as an 
agent for initiating a turbulent boundary layer. There 
is further the important problem of buffeting caused by 
the long bubble type of separation; although this has 
only just begun to be studied the results are already 
finding practical application. 

The effect of compressibility on long bubbles and the 
question of transonic reattachment are also important, 
Also Chap 


man ef a/.*? have made a preliminary general study of 


and have been touched on in reference 51. 
separated flows in supersonic (and subsonic) streams 
The work of Holder and Gadd 
significant contributions to these problems, in which it 


and of Pearcey*' are 


may be expected that the general approach of th 
present paper will be helpful 

Finally, a hypothesis concerning the mechanism of 
the leading-edge stall has emerged from the physical 


rhis 


helps to explain the scale effect on the stalling charac 


picture of laminar separation bubbles given here 
teristics of airfoils experiencing bubble separations, but 
it is still not possible to calculate the maximum Iift 


achieved in a particular case 


) aia 
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By using the principle of least work, it is shown that it is possible 
determine the shape of the cross section of the tube as a fune :' 
‘| tion of longitudinal curvature and als ifter the shape is A1NTRODUCTION 
. of found, t letermine the bending moment carried by the tubs , 
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ber With this latter available, bending modulus of rupture charts for : ; 
und tubing are constructed which appear to agree well with deform during bending and that this deformation ts 
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Fic. 2 [ypical stress-strain curve used in this investigation. 


stressed above the elastic limit or in the plastic range 
for the material. 

During ovalization, it is assumed here that tubes 
which are initially round become elliptical in shape, that 
this circumferential deformation is inextensional, and 
that lines originally normal to the middle surface remain 
so after loading; these are the same assumptions used by 
Chwalla. 


subjected to pure bending and is of sufficient length so 


Other assumptions are (1) that the tube is 


that the end conditions have a negligible effect on the 
strength at its center and (2) that the material is iso 
tropic and has the same stress-strain curve in tension 


and compression. 
ANALYSIS 


Work of Deformation 


In order to establish a theory, a way must be evolved 


to determine the amount of ovalization of the cross 


section of the tube for a given longitudinal radius of 
curvature. It will be shown that this can be done by 
using the principle of least work or, in other words, by 
minimizing the total internal work of the bent tube. 
In order to show that the method used to determine the 
internal work applies in both the elastic and _ plastic 
ranges, an expression for the former will be derived 
we see that in the elastic 


first. Referring to Fig. 1, 


range the work per unit length is 


dy (1 2)orpe, + (1/2) orer] dé ds 


where o,, «,, or, and =p are, respectively, the longi 
tudinal and transverse stress and strain. 
Substituting for o, and o, their values as functions 


of £, M, and € 


daly’ [A /2(1 — wu?) | [Cen + werde, + 
(er + pez)er |dtds (1) 
= [EF 2(1 — wp?) |[en? + er? + Quezer |déds 
This may be written 
dW = (E/2)(e;?)déds (2) 


T 


C 
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Ii 
e ¥Y I1T/Q u-)i(exn° + er* + Queze 
Eq. (2) may also be written 
dy (1 2)(Ke,edéds | 2)o,e,d¢ds 
where go; and ¢, are coordinates of the stress-strain cury, 


for the material. Since (1 2)¢,e,; equals the area unde 
the stress-strain curve, it is readily seen that the worl 
per unit volume at any point in a two-dimensional 
stress system can be determined by first computing ¢ 
using Eq. (3) and then determining the corresponding 
area under the stress-strain curve for the material. § 
far this has only been proved to hold for stresses stil] 
within the elastic range for the material. 

In the plastic range let us begin in the usual manner 
2. then reduces to th 


by assuming yu l Eq. (3 


following: 


which is the equation for the strain intensity in the 
the 
If it can be shown that the corresponding 


plastic range in accordance with deformatior 


theory. 
area under the stress-strain curve also gives the work 
per unit volume, then the theory will be continuous i1 
the elastic and plastic ranges. This is important her 
since a section can at the same time have stresses 1 


both ranges. Let us begin by assuming that this is so 


or that (see Fig. 2 
1” JS 0 de JS Eye de 
J I,(2 V3)(en> + €7 t e,er)'/? X 
((1/YW3)(eLr> + er? + exer x 
(Verde T end €7 e, dey T erde ; 
(4.5 JS Es le + (1/2)e7 de T 
ep + (1/2)« der; 
S (a de ord ey 
since” a (4 3)F, len + (1/2)e7 
and O71 (4/3 E | €7 1/2)e 
Eq. (5) represents the total work per unit volume 
This means that the work per unit volume can be 


determined in both the elastic and plastic ranges by 
taking it equal to the area under the stress-strain curve 
This 


has another advantage in that a varying value of yu can 


corresponding to the value of ¢, given by Eq. (3 


be used which can be important in elastic-plastic cases 
since w for some materials might vary by as much as 
100 (0.25 to 0.50 from 
elastic to plastic range). The total work may then be 
determined by computing the work per unit volume at 


per cent as stress increases 


various points using the above process and integrating 
over the total volume of the material 


Application of Principle of Least Work 


In what follows it will be shown that this principle 
applies in both the elastic and plastic ranges for the 


problem being considered. Let us take a round tube 


in 
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ds 


rain ¢ urve 
rea under 
the work 
nensional 
Ilputing ¢ 
‘sponding 
arial. § 


esses still 


] manner 
‘Ss to the 


v in the 
TInation 
sponding 
he work 
NuUOUS In 
ant here 
resses in 


11S 1S SO 


‘olume 
can be 
ges by 
1 curve 

This 
f uw can 
C Cases 
uch as 
; from 
hen be 
ime at 
rating 


nciple 
yr the 
| tube 


BENDING STRENGTH OF 
The 


and bend it into a fixed circular arc of radius p. 
stable 


tube will ovalize and assume a position of 
equilibrium. Since the assumed cross section is ellip- 
tical (see above), let us further assume that the semi- 
major axis varies by a small amount. With the ends 
fixed against rotation no change in the external work 
occurs. The internal work therefore remains the same, 


and we can write the following equations: 
JSAWdV = Af WdV = 0 


is the work per unit volume and | 
If the shape of the ellipse is de- 


where IV is the 
volume of the tube. 
fined by the quantity a 7, the ratio of its semimajor axis 
and the radius of the original round, the above equations 
signify that at this point of stable equilibrium the 
change in internal work due to a small change in a/r is 


equal to zero, or that 
OW O(a r) = O 


[his means then that at the stationary point the total 
internal work with p held constant and a r varied must 
be either a maximum or minimum. Reference to IV vs. 
a ry curves which were obtained by a process of numer- 
ical integration (see Fig. 4) indicates that this point is 
always a minimum. More extensive plots of IV vs. ar 
also show that no other maximum or minimum points 
for which the slope is zero occur, at least within the 
physical limitations of the problem. 

Hill” extremum and 
principles for the plastic-elastic range which could be 
These would indicate 


also discusses variational 
applied to a certain extent here. 
of course that the absolute minimum of total internal 
work would be realized by use of the proper cross- 
sectional shape. The assumption of an ellipse is, how 
ever, a convenient one and apparently sufficiently 


accurate for the purposes here. 


Application to Round Tubing 


This general theory will be applied to determine the 
bending strength of round tubing. Since the stresses 
in all transverse sections are identical, we only need to 
consider a tube of unit length; and also, because of the 
section, the 


symmetrical condition of a_ transverse 


integration need be performed only over one half of it 


Y 
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Quadrant of ellipse showing sections and points used in 
numerical integration 
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Fic. 4. Curve illustrating method for obtaining a + correspond 


ing to given value of p using principle of least work 


The 


section is not quite symmetrical about the neutral axis 


(one quarter above and below the neutral axis) 


owing to small differences in the states of stress and 
strain, although the resultant forces above and below 
the neutral axis are practically identical for the case 
examined. Even with these simplifications, however, 
no closed solution for the required integration was 
found, and therefore a numerical integration process 
was resorted to. In this, the one-quarter circumference 
of the ellipse was divided into six equal parts (1-7 
Elliptic integrals were used in order to secure 
sufficient accuracy. At each of these division points 
the wall was divided into four equal parts. Taking 
section | for example, the work per unit volume was 


Fig. 3). 


computed at each of the five division points, and the 
integration across the wall was performed by use of 
Simpson’s rule. After this process was completed at 
each of the remaining sections 2-7, the circumferential 
integration was carried out, again by use of Simpson's 
rule. This process of course was repeated for a quad 
rant below the neutral axis and the total work found by 


adding the results and multiplying by two. 
Deformation Relationships 

The data to be substituted in Eq. (3) are as follows: 
P)\Vn T+ ¢ COS @ (()) 


e, = Vp = (1 


where p is the longitudinal radius of curvature of the 
bent beam, and a is the slope angle of the middle sur 


face. 
l 
cos a = t : 
1 -+ (dy dx)*|!- 
l ss 
& 3 aero , 
»1 + (67/a?)[(1 — ¥m?/b?)/(9 h*) |; 
} 
1 CCL. pr) (i/f}i — pity 
ClU/pr) — (r)] — uw p ‘) 


Here pr and r are the mean radii of curvature of the 


ellipse and circle, respectively, and €; 5 p is the 
mean or middle surface strain for the ellipse 
ly pr = (d*y/dx*)/{1 + (dy/dx)?]*/? = 
ab/[b? — vm? + (a? b*)y i) 
a=) (1 + 6)/(11 +e) Jt 10 
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Fic. 5. Curve illustrating method of obtaining maximum 
moment 

where 6 = (1 + ep,)(1 — mwprepz)*? — 1 

For Ci < €pl, b= Up! 


This value of Poisson's ratio u appears to check well 
with experimental results.° Since the values of uw and e; 
are interdependent, the latter have to be determined 


by trial. 


Determination of Ovalization 


The amount of tube ovalization corresponding to a 
given value of longitudinal curvature p also has to be 
determined by trial and is accomplished as follows: 

(1) Assume a value for p. 

(2) Assume a certain amount of ovalization or, in 
other words, a value for ar. 

(3) With this and other information given above, the 
value of ¢; at any point may then be determined; and 
the corresponding work per unit volume may be deter- 
mined by use of the following equations for the stress 


strain curve 


(11) 


E){1 + (3/7)(6;/ 00.7)" ~'] 


1) ](o;/o0.7)"~'; (12) 
Eq. (12) is the area under the stress-strain curve and is 
obtained from Eq. (11) by integration. 

(4) The total work for the deformed section is ob- 
tained by the numerical integration process as previ- 
ously explained. 

(5) Since the correct value of ovalization for a given p 
is obtained when the work is minimum, this process is 
repeated until the minimum is obtained. In doing this 
p is held constant and a r varied, and the results are 
plotted as a curve of Il’ vs. a r as shown in Fig. 4. 

Obviously the above process is too laborious for hand 
that it 
gramed and accomplished using a fairly high-speed, 
The data shown 
for example, was computed using a 


calculation, but it was found could be pro- 


high-capacity digital computer. 
plotted in Fig. f, 


Sendix G-15 in a matter of approximately | hour. 


Transverse Deformation Relationship—Tubes 
Other Than Round 


This method can also be used for cylindrical sections 
which are not initially round by substitution of the 
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following equation for er: 
er = [(1/er) — (1, pr’)] — w(m/p) 


Here (1. pr) (1 pr’) represents the change in trans 


verse curvature. This equation can also be used for 
tubes which are supposed to be initially round but 
which have some ovality. For convenience, the initial 


cross-sectional shape would again be assumed elliptical 


Determination of Bending Moment 


Once the ovalization corresponding to a given radius 
of curvature is determined, the next step is to deter 
This 
is done by summing up the moments of the longitudinal 


mine the moment carried by the ovalized tube 
stresses about the neutral axis. The equation for the 


longitudinal stress is as follows: 


(1 - a”) |(e + mer) 15 


where uw, €,, and ey are as given by Eqs. (6)—(10 
iD = a;/e 
strain curve for the material. ¢, is determined by Eq 
After the dM/ 


is obtained at each point (see Fig. 3), the total moment 


and can be determined from the stress- 


(3) as previously explained. o,yvd¢ds 
for the section is obtained by the same numerical 
integration process used for the determination of the 
total work. 

In order to determine the maximum moment that a 
given section can carry, several values of moment must 
be determined for varying values of p. As the tube 
bends and p decreases, at some point the moment 
carrying capacity of the tube reaches a maximum and 
then starts to decrease due to the excessive ovalization. 


A curve illustrating this is shown in Fig. 5. 


DISCUSSION OF RESULTS 


Ultimate Bending Strength of Round Tubing 


A study of this subject indicates that the mode o 
failure varies and that it can be divided into ranges in 
accordance with the diameter to thickness or D ¢ ratio 
somewhat as follows: 

(1) AtaD/t = 
failure with some, but not appreciable, ovalization or 
The strength can 


2 (solid round) the failure is a material 


other distortion of the cross section. 
be computed conservatively by the method of reference 
ff 

(2) In the range of D/t's from 3-10, the failure is still a 
material failure but takes place with some ovalization 
of the cross section. The amount of ovalization in 
creases as the D ¢ increases. 

(3) In the range from 10-50, failure is by excessive 
ovalization in the manner already explained. It should 
be noted that in a test the apparent mode of failure in 
this range would be either material failure or failure by 
buckling, since one or the other of these would certainly 
occur after the maximum moment had been reached 
Some carefully run tests are needed in order to deter 


mine this point accurately. 
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BENDING STRE 


1) In the range of D ?¢’s > 50, failure may be by 
buckling which takes place after some ovalization of the 
eross section. In order to indicate the buckling type of 
failure, however, rather large imperfection factors* 

have to be that the initial 


imperfections have quite an effect on the failure mode 


used, and this indicates 


ind strength of the tube. A more accurate theory for 
the elastic and plastic buckling of bent and ovalized 
tubes is needed in order to secure a better check on this 
point \ method of computing the buckling curve 
will not be given here since it is quite an extensive 


itself and one which needs considerably 


problem in] 
more investigation and refinement. 

Fig. 6 shows curves for 260,000 ultimate tensile 
strength steel which were plotted in accordance with 
items (1)-(3) given just above. In accordance with 
item (3) then, for those sections which do not fail by 
material failure, excessive ovalization rather than 
buckling is considered to be the cause of failure. The 


the stress-strain curve for 
29 000,000, 


plotted is the com 


constants used in Eq. (11), 
the material, were as follows E m = 
16.0, and o 214,000. The F, 
puted ultimate moment divided by the section modulus 
The failure with no 
the 
ovalization 


of the original round cross section. 
ovalization curve was computed by method of 
reference 7. The 
determined by computing data for tubes with D/t 
30, and 50. 


failure with curve was 


values of 10, 


Bending Modulus of Rupture Curves for Round Tubing 


The two curves shown in Fig. 6 will form a smooth 
curve since the ovalization is gradually and continu 
ously the 
general, gradually decreasing as the D ¢ increases. By 


increasing and longitudinal stresses, in 
taking the curves shown in Fig. 6 and plotting them as a 
continuous curve (see Fig. 7), there is obtained a curve 


which can be used to obtain the bending modulus of 
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Computed and experimental bending modulus of rupture 
data—260,000 U.T.S. round tubing 


rupture for the material in question. Similar curves 
for a number of different steels have been plotted, and 
the agreement with available test data has in all cases 


been about the same. 


Stress-Strain Curves 


In plotting these curves, the conventional rather than 
the so-called true" stress-strain curve was used. Since 
the latter is higher, particularly at the ultimate strength 
end, the corresponding bending modulus curve would be 
higher at the left end but about the same at the right 
end. This would tend then to rotate the curve so that 
it would appear to lie closer to a mean curve for all the 
test data. This was not investigated, however, for the 
following reason: 
to be used as design curves, and many designers or stress 


the curves were being constructed 


analysts would probably not use them if they found that 
some of the available test data plotted below the curves 


CONCLUSIONS 


(1) The assumptions and method of numerical in 
tegration employed apparently give sufficiently accu 
rate results for practical purposes. 

(2) Use of the deformation theory again appears to 
lead to sufficiently accurate results. 

) In some test data, a failure moment smaller than 


\- 


the maximum moment taken during the test might 
readily be reported. This is due to the fact that, 
because of the ovalization, the moment can reach a 


maximum and then start to decrease before the more 
type of failure, by fracture or buckling, 
occurs. Also, test data 
moments than those obtained by the theory due to the 


apparent 
some will indicate higher 
effect of some local restraint. 

(4) The use of the true stress-strain curve would 
probably give a bending modulus curve more in line 
with the mean test curve, but it would not be generally 
acceptable for design use because some of the test 


points would fall below the curve. 
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(5) The methods used to construct the bending 
modulus curves appear to give results in line with the 
available test data even for some of the materials which 
are known to be somewhat anisotropic e.g., 260,000 
U.T.S. steel. 

In order to check more completely the methods 
employed and conclusions reached, there is need for a 
more accurate theory for the elastic and plastic buckling 
of bent and ovalized tubes and also for several carefully 
conducted tests run with some of the above-mentioned 


points in mind. 
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Aerodynamic Action of Triangular Horn-Balanced Control Surfaces . 


(Continued from page 578) 


considered by Coale enables the final formulas to be 
reduced to functions of just one parameter, the reduced 
aspect ratio of the balance. This simplification could 
not be achieved in the case of a triangular balance. 
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A Simple Model Study of T 


— 


Temperature and Thermal Stres 
Distribution Due to Aerodynamic He ating’ 


GABRIE 


SUMMARY 


The present work is concerned with the determination of tran 
jient temperatures and thermal stresses in simple models intended 
simulate parts or the whole of an aircraft structure of the built 
ip variety subjected to aerodynamic heating 

The first case considered is that of convective heat transfer 
ito one side of a flat plate, representing a thick skin, and the 


fect of the resulting temperature distribution in inducing ther 


mal stresses associated with bending restraint at the plate edges 
Numerical results are presented for the transient temperature dif 
in the plate when the environment temperature first 


the 


entials 


linearly with time and then remains constant, 


mcreases 


period of linear increase representing the time of acceleration of 


he aircraft Corresponding thermal stress information is pre 


sented 


rhe second case is that of the wide-flanged I-beam with convec 


tive heat transfer into the outer faces of the flanges. Numerical 


results are presented for transient temperature differentials for a 


wide range of values of the applicable parameters and for an en 


nment temperature variation as described above Corre 






sponding thermal stresses in a beam of infinite length are deter 
nined 


I finite 


\ theoretical analysis of the stress distribution in a beam 


length is carried out and numerical results obtained for 


case \n experimental investigation of temperatures and 


stresses in such a beam is described, and results are presented 


which indicate good agreement with corresponding theoretical 
results 
SYMBOLS 
semiwidth of I-beam 
( specific heat of plate or beam material 
f(g nondimensional temperature distribution in 
plate, (JT — 7Ty=0)/A7 
/ heat-transfer coefficient at heated surface 
h heat-transfer coefficient at joint between beam 
flange and web 
R thermal conductivity of plate or beam material 
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AAA 


> 


semilength of I-beam 

semidepth of I-beam 

time 

time of linear increase of environment tempera 
ture 

displacement in x-direction 


coordinate in plate representing distance from 


insulated surface, positive toward heated 
surface 

coordinates of I-beam (Fig. 8 

cross-sectional area of flanges and web. re 


of I-beam 


coefficients in series for 


spectively, 

\iry stress function for 
flange 

coefficients in series for Airy stress function for 
web 

Young’s modulus of elasticity 

stress correction factor for 


compression (Eq. 5 


stress correction factor for tension (Eq. 6 

modified stress correction factor for 

Eq. 7 

modified 
Eq. 8 


compres 
S1on 
stress correction factor for tension 
temperature in plate or beam 


temperature at middle of beam web 


differential between surfaces of 


temperature differential between mid 


temperature 
plate, 
dle of web and edge of flange 

initial temperature in plate or beam 

environment temperature 

reference environment temperature 

strain energy 

coefficient of thermal expansion 

ratio of flange thickness to web thickness, 6,/6 

plate thickness 

nondimensional beam coordinate, 2 

nondimensional beam coordinate, y 

thermal diffusivity of plate or beam material, 
k/ pC. 


} 


nondimensional heating parameter hé/k in 
plate, hs?/ké,; in I-beam 

nondimensional heat-transfer coefficient at 
web-flange joint, A;s/k 

2+ Pp 

Poisson's ratio 

nondimensional beam coordinate, x/6 in plate 
x/s in I-beam in Section (4), x// in I-beam 


in Sections (5) and (6 


density of plate or beam material 
normal stress components 


intensity of normal force applied on end of 


beam 


nondimensional time variable, «t/é? in plate, 


xt/s? in beam 
of linear of en 


nondimensional time increase 


vironment temperature 








612 JOURNAL OF THE AERONAU 
T9 
r — 
tr at (Ad) maz, —9 
Ts Tyas Tx = shear stress components 
i) = nondimensional temperature (7,, — 7 
T., — T,), Airy stress function 
do, = nondimensional environment temperature, 
(Pa = Ee T., -— 1 
Ao = AT/(T., — 1 
y = ratio of I-beam length to width, //> 
w ratio of I-beam width to depth, 6/s 
Subscripts 
flange 
1 = web 
(1) INTRODUCTION 


N THE RAPID ACCELERATION of an aircraft to very high 
I speed, friction in the boundary layer may raise the 
temperature of the surface of the aircraft with sufficient 
rapidity that large temperature differentials may oc- 
cur in the structure before interior portions heat up. 
This may result in thermal stresses of considerable 
magnitude. 

The heat transfer into the surface of the aircraft may 
be treated in a manner similar to convective heat trans- 
An appropriate heat-transfer 
coefficient with the difference 


between the temperature at the solid surface and the 


fer in low speed flow. 
is used in conjunction 
so-called ‘‘adiabatic wall temperature,’’ or temperature 
which would be attained in the air at the solid surface 
if that surface were insulated.! The heat-transfer 
coefficient and adiabatic wall temperature depend upon 
the properties of the ambient air, the velocity of flow 
outside the boundary layer, and the thickness and type 
of boundary layer (laminar or turbulent). 
Determination of the temperature 
thermal stress distribution in the case of a typical built- 


transient and 
up aircraft structure is exceedingly difficult and cum- 
bersome. Much of a 
however, from the analysis of simple structures repre- 


general nature can be learned, 


senting parts of the built-up structure or constituting 
simple models intended to represent the whole struc 
ture. 

Thus, in the present work, attention is concentrated 
on the flat plate with convective heat transfer into one 
face and the wide-flanged I-beam with convective heat 
transfer into the outer faces of the flanges. The former 
is intended to represent the thick skin of a high-speed 
while the latter is intended to represent a 
The interior 
surfaces of the structure are assumed to be insulated. 


aircraft, 
simple type of built-up wing structure. 


Heat transfer by radiation or convection between parts 
of the structure is neglected. 

Studies of this type have been reported in a number 
of publications in recent years,° particularly on the 
of infinite 
have been presented for selected values of the 
In the present work, 


wide-flanged I-beam length, and results 
signifi- 
results 
are the 


parameters reduced to nondimensional form, and the 


cant parameters. these 


extended to cover large ranges of values of 


problem of thermal stresses in a wide-flanged I-beam of 
finite length is considered in some detail. 


T 


rc 
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(2) TRANSIENT TEMPERATURE DISTRIBUTION IN 4 
FLAT PLATE 
The case of a flat plate insulated on one side with 


convective heat transfer on the other side is first con 


sidered. The following assumptions are made 


(1) 
coefficient are uniform over the heated surface 


The environment temperature and heat-transfe 
the lat 
ter being constant with time. 


(2) The temperature in the plate is initially unj 
form. 
(3) There is no heat transfer through the edges oj 


the plate, so that the flow of heat takes place only in a 
direction normal to the surface of the plate. 

(4) Variation of the thermal properties of the ma 
terial of the plate with temperature may be neglected 

An equivalent problem for the special case of con 
stant 7; 
and its solution is well known. 
infinite plate at initial uniform temperature 7 
to at 
on both sides. Heat flows into the sur 


is a familiar one in the field of heat transfer 
It is the problem of an 
, sub 
jected suddenly an environment constant tem 
perature 7; 
faces of the plate by convective heat transfer. 

This problem is seen to be identical with the present 
problem if, in the former, only half the thickness of the 
because of there 


plate is considered, since, symmetry, 


is no heat flow across the middle surface, which maj 
thus be considered as an insulated surface. 

The classical solution to this problem is given as an 
infinite transcendential series.”,» Numerous charts and 
tables presenting results for various values of the param 
eter \ are available in the literature 
115, 116, 


These results are used here to plot a curve of 


(see reference 2 


and reference 3, p. 
(Ad) maz 


that is, the nondimensionalized maximum tempera- 


129 


pp. 274, 285, and 


ture differential across the two faces of the plate during 
the heating period the The 
nondimensional time 7 at which this maximum occurs | 


versus parameter A. 


has also been plotted. These plots are shown in Fig 


1. It will be shown in the next section that this maxi 


mum temperature differential provides a measure of the 
maximum thermal stress in the plate. 

The assumption of constant 7), and h is not realistic 
for the case of aerodynamic heating, since it implies an 


infinite acceleration, that is, an instantaneous achieve 
ment of high-speed flight. In reality, the period oi 
acceleration may be of the order of several minutes 
during which time the adiabatic wall temperature 
represented here by the environment temperature, and 
the heat-transfer coefficient are increasing. The m 


troduction of a variable heat-transfer coefficient would 
complicate the problem very greatly, however. On the 
other hand, the introduction of a variable environment 
temperature does not complicate the problem unduly ! 
Re 


cause of the linearity of the governing differential equa 


its variation is assumed to be of a simple form 


tion, such a variation of 7, may be taken into account 
by the application of Duhamel’s integral in conjunc 
tion with the time history of temperature in the plate 
following a sudden change in environment tempera 


pe 


th 
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. \ y . . . y "yy 
an Cg giving a compressive stress ak AT7/2(1 v) at the 


side with 


S first con 


ide 


at-transfer 


ce, the lat 


itially uni 


le edges oj 


e only ina 


~ the ma- 


neglected 
se of con 


t transfer 


blem of an 


e T';, sub 
tant tem 
0 the sur 


he present 
less of the 
try, there 
hich may 


ven as an 
harts and 
he param- 
ference 2 
p. 122 
ft (Ad mer 
tempera- 
ite during 
A. The 
im occurs 
m in Fig 
his maxi 
ure of the 


t realistic 
mplies an 
achieve 
period of 
minutes 
perature 
ture, and 
The in 
nt would 

On the 
ronment 
anduly 1! 
‘m. Be 
ial equa 
account 
conjune- 
he plate 
empera- 


| 


ture. While the problem may be handled analytically 
in this manner, it was decided, after considerable inves- 


tigation, that the use of a finite difference numerical 


method 11] 
computing machine would be more expeditious in the 
Consequently, 


conjunction with an automatic digital 
generation of a large volume of data. 
. finite difference procedure was set up. 

The time history of 7, was assumed to have the form 
shown in Fig. 2 and is represented mathematically as 
follows 


for ON t< £8 


7 fr, ta Tt, 
a for 


t> t, { 


A program of computations was carried out on the 
Whirlwind I automatic digital computer at the Massa- 
chusetts Institute of Technology for a variety of values 
of the parameters \ and 7,. The effect of these param- 
eters on the maximum temperature differential at- 
tained in the plate is shown in Fig. 3. 

In addition to the maximum temperature differen- 
tial, it is also important, for the purpose of thermal 
stress determination, to know the distribution of tem- 
perature through the thickness of the plate. This 
distribution may be expressed conveniently in terms of 


the nondimensional function, 
f(t) = (T — T,.9)/AT (2) 


The computed results indicate that, at the time of 
maximum temperature differential, the function /(£) 
approaches the parabolic form, f(£) = &, for low values 
of \ when 7, varies as a step-function, and for all 
values of \ when 7°, has an initial linear gradient and 
the nondimensional gradient time 7, is sufficiently 
It is thus convenient to adopt this parabolic 
large 


large. 


distribution as a reference distribution. At 
values of \ and zero or small values of 7,, the tempera- 
ture gradient at the heated surface is more pronounced 


than that given by the parabolic distribution. 


(3) THERMAL STRESSES IN A FLAT PLATE 


The thermal stresses resulting from transient tem- 
perature distributions in a flat plate, as discussed in 
the preceding section, are now considered. 

The usual assumption in plate theory—namely, that 
stresses in planes normal to the plate may be neglected 

is made here. Furthermore, it is assumed that the 
plate is free to stretch or contract in its middle surface, 
but is restrained against rotation at its edges, and that 
a and / do not vary with temperature. Under these 
circumstances it is easily shown that the thermal stres 


ses in the plate are given by 


=o lakAT/(1 — v)] X 
1 
' stat 
y JAS) — | E)dé | ) 


For a linear temperature distribution, Eq. (3) be- 


comes 


—[aKkAT/(1 — v)|[é —a 2)] (4) 


o oO = 


heated surface and an equal tensile stress at the insu 
lated surface. 

Such a temperature distribution is not consistent 
with the boundary condition for heat flow at the in- 
sulated surface, and, as seen in the preceding section, 
the actual temperature distribution for the type of 
heating considered will be approximately parabolic 
in most cases, and the stress will be as shown in Fig. 4. 
It is, however, convenient to use the linear distribution 
for the purpose of computing a reference thermal stress 
to which a correction factor may then be applied for 
other temperature distributions. Thus, stress correc 
tion factors for compression and tension, respectively 


are defined as follows 


{ ’ 
K Z l f(E\dé ) 
' Jo y 
*! 
Ks 2 f(ejdé 2 K 6 
/ 0 
For a parabolic distribution of temperature |/(¢ 
K, = 4/3, K 2/3 


When the temperature gradient at the heated sur 
face is still more pronounced, A, will be larger, with 2 
as its limiting value, while A, will reduce toward zero 

The transient temperature data obtained from the 
digital computation program have been used to obtain 
transient stress data. Fig. 5 shows the variation ot 
K, and K, with time for a typical case. 
this variation, the maximum stresses in the plate are 


Because of 


not reached simultaneously with the maximum tempera 
ture differential. 

In order to express the maximum stresses in terms ot 
the maximum temperature differential, it is desirable 


to define modified stress correction factors as follows 
Ke’ = 1A6/(Ad)ece| Kz (7 


K,’ = [Ad/(A6),.-1 Be (S 
The maximum values of these modified factors, when 
multiplied by the maximum temperature differentials, 
will give the maximum stresses. Their variation with 
time is shown in Fig. 5. It is seen that the maximum 
compressive stress is reached somewhat earlier and the 
maximum tensile stress somewhat later than the maxi 
mum temperature differential. The maximum value 
of K,’ is 1.55, and consequently considerably greater 
than the value (1.333) of A, for a parabolic tempera 
ture distribution. On the other hand, the maximum 
value of K,’ is 0.64 and, consequently, only slightly less 
than the value (0.667) of AK, for a parabolic distribu 
tion. 

In general, for large values of the parameter \ and 
r, = 0, the maximum value of K,’ tends to be substan 
tially larger than the value of A, for a parabolic tem 
perature distribution, while the maximum value of 
K,’ tends to be slightly less than the value of A, for a 


parabolic temperature distribution. As A decreases, 
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SIMPLE MODEL STUDY OF 


the times at which KA,’ and A,’ reach their peak values 
appr ach the time for (A@) maz, and the peak values of 
these correction factors approach the corresponding 
values for a parabolic temperature distribution 

namely, 1.333 and 0.667, respectively. These results 
are shown for compressive stress in Fig. 6, along with 
the variation of K, at (A@) ma; with X. 
tensile stress, the maximum value of A, is either equal 
to or slightly less than 0.667 throughout the range con- 
sidered, and may thus be taken as 0.667 in any practical 


In the case of 


computation. 

The effect of 7, on the value of AK, at (A@) mar is Shown 
in Fig. 7. It is seen that A, at (A@)naz approaches the 
value 1.333 as 7, increases. The effect of 7, on the 
maximum value of K,’ was not obtained. It is reason- 
able to expect, however, that this quantity will ap- 
proach the value 1.333 as 7, increases. The maximum 
value of A,’ can be expected, in all cases, to remain 


close to the value 0.667. 


1) TRANSIENT TEMPERATURE DISTRIBUTION IN A 
WIDE-FLANGED I-BEAM 


A simple model of a built-up type of wing structure is 
the I-beam with wide flanges. The present section is 
concerned with the transient temperature distribution 
associated with the heating of such a beam by convec- 
tive heat transfer into the outer surfaces of the flanges. 

The following assumptions are made. 

|) The environment temperature and convective 
heat transfer coefficient are uniform over the outer 
surfaces of the flanges, the latter being constant with 
time as well. 

(2) The temperature in the beam is initially uniform. 

(3) The flow of heat takes place entirely in directions 
parallel to planes representing cross sections of the 
beam. That is, the heat loss out of the ends of the 
beam may be neglected, so that the temperature dis- 
tribution at any instant is the same in all cross sections 
of the beam. 

(4) The flanges and web are sufficiently thin that 
temperature variation through their thickness may be 
neglected. 

(5) The inner surfaces of the flanges and both sur- 
faces of the web are insulated, so that heat transfer be- 
tween parts of the beam takes place by conduction 
only. 

(6) There is no heat transferred through the edges of 
the flanges. 

(7) Variation of the thermal properties of the ma- 
terial with temperature may be neglected. 

The problem thus reduces to one of one-dimensional 
heat flow in the flange and in the web. 

A simplified treatment of this problem has been given 
by Hoff and Torda,' in which the heat loss from the 
flanges into the web is neglected in the determination 
of the transient temperature variation in the flanges, 
and this temperature variation is then applied as a 
boundary condition on the web. More refined analyses, 
based on assumptions which are the same as those 
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made here, have been given by Pohle and Oliver,' 
Schuh,* and Parkes.’ In the present work, the results 
of this previous work are extended considerably and 
information is provided on the effect of a finite period 
of build-up of the environment temperature 

The present analysis satisfies the transition conditions 
at the joints between the flanges and web exactly 
The beam is assumed initially to be built-up, in the 
sense that the flanges and web are not integral but are 
connected by screws or pins or some similar arrange 
ment, leaving a surface between the flange and web 
across which there may be effectively a temperature 
discontinuity. The solution for a beam with integral 
flange and web is then obtained as a special case. 

The coordinates and dimensions of the beam cross 
section are shown in Fig. 8. The depth of the beam is 
measured between the middle surfaces of the flanges. 

The partial differential equations in nondimensional 
form for the temperature in the flanges and web are, 
respectively, 

Og, Or = (0°¢,/O0E*) — Ald — G, 8) 
0¢,/OT = O°hy/ On” (10 

Taking symmetry into account, it is necessary to 

consider the flange only in the range 0 < & < b/s, and 


the web only in the range 0 < n < 1. The boundary 
and initial conditions may now be written 


(1) até = b/s, O¢g,/0E = O 
(2) aty = 0, O¢,/On = O 
(3) até = 0,7 = 1, 


27(0¢,/OE) = Aj(Gr — Hw) Od, On 


(4) atr=0, @ = d& = 1 


As in the case of the plate, a finite difference numeri- 
cal method was used in the solution of Eqs. (9) and 
(10). 

A program of computations was carried out on the 
Whirlwind I computer for a variety of values of the 
parameters A, y, and 7,, the time history of the environ- 
ment temperature being assumed to have the form 
shown in Fig. 2. The influence of the parameters 
\ and y on the maximum temperature differential in 
the beam is shown in Fig. 9 for the case of 7, = 0. 
The width of the flange was taken to be twice the depth 
of the web in all cases. The web and flange were as- 
sumed to be integral. 

The influence of the gradient time, 7,, is shown in 
Fig. 10 for the case of y = 1. The results are expressed 
as a factor modifying the results obtained for the case of 
tr, = 0. The effect of 7, is seen to be greatest at inter- 


mediate values of X. 


(5) THERMAL STRESSES IN A WIDE-FLANGED I-BEAM 
OF INFINITE LENGTH 
The thermal stresses due to the transient tempera- 
ture distribution in the I-beam described in the pre- 
ceding section are now considered. 
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The beam is first assumed to be infinitely long. Since 
the temperature distribution is the same in all cross 
sections of the beam and is symmetrical about the x and 
y axes (see Fig. 8), it follows that cross sections must 
remain plane and parallel during heating of the beam. 
The assumption of thinness of the flange and web and 
the fact that the temperature is uniform through their 
thickness results in a condition of plane stress in each. 
As the beam is free to deform within any cross sec- 
tion, the only stresses present will be in the longitudinal 
direction, that is, ¢, = ¢, = 0. It is easily shown® that 
under these circumstances the longitudinal stress is 


given by 


oz, ak(AT) maz [1 (Ad) mar |) pl E) — Aj (11) 
o, /a@E(AT) mar = [1/(Ab) mar liGu(n) — A} (12) 
where 
ye 1 
A : | o(édé + 
<Yw T l J 0 
l I 


du(n)dn (15) 
2yw + /J0 


It should be noted that & is defined differently here 


than in the preceding section. 
It is seen that 


ak (AT) maz (14) 


LCs Janu — lz Jiud 


The temperature data of the preceding section have 
been used to compute stress data for the wide-flanged 


I-beam. In addition to the parameters \, y, and 7,, 
the parameter w is now significant. While the tempera- 
ture distributions were computed for w = 2, it is in 


most cases sufficiently accurate to use these distribu- 
tions for other flange widths, cutting the distribution 
curve off at the appropriate point for a narrower flange 
and extending it horizontally for a wider flange. This 
procedure introduces the greatest error at low values of 
\, where the effect of heat loss to the web is felt far out 
on the flange. 

Fig. 11 shows, for 7, = 0, how the stress range in the 
beam at the time of maximum temperature differential 
divides between maximum tensile stress and maximum 
compressive stress. The maximum tensile stress oc- 
curs at the middle of the web, while the maximum com- 


pressive stress occurs at the edge of the flange. The 
curves are identified in groups of constant ratio of flange 
to web cross-sectional areas, A;/A,, = 2yw. 


The effect of gradient time on the maximum values 
of ¢,/ak(AT) naz was found to be slight in all cases. In 
the extreme case of } = 100, the maximum tensile 
stress is increased about 1.5 per cent and the maximum 
compressive stress is decreased about 8 per cent when 
the relative gradient time (7,), is increased from 0 to 8. 
In any practical application the stress may be deter- 
mined by the use of Fig. 11 in conjunction with the 
temperature data of the preceding section. 
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(6) THERMAL SPrRESSES IN A WIDE-FLANGED I-Beay 
OF FINITE LENGTH 


The effect of free ends on the thermal stress distriby 
tion in the wide-flanged I-beam is now considered 
The assumptions regarding temperature distributior 
remain unchanged. The coordinate axes of the beam 
are shown in Fig. 8. The origin is at the center of th 
beam. 

The flanges and web are again assumed to be suf 
ficiently thin that stresses normal to their surface may 
be neglected. Stresses in the plane of the flange and 
web are considered to be uniform through their thick 
nesses. Deflections of the flange and web normal to 
their respective planes are assumed not to occur. This 
precludes buckling of the flanges under compressive 
stress. 

On the basis of these assumptions, the problem r 
duces to one of plane stress in the flange and web, with 
the satisfaction of boundary conditions at the joints 
connecting them and at their free edges. 

The stress distribution may be found conveniently 
as the superposition of two component stress distribu 
tions. The first component stress distribution is that 
corresponding to the beam of infinite length. The 
second component stress distribution is that necessary 
to result in satisfaction of the boundary conditions on 
the ends of the beam, and is obtained by applying to 
the ends of the beam a longitudinal load distribution 
equal and opposite in sign to the stress distribution in a 
cross section of the infinite beam. Since this load dis 
tribution is self-equilibrating, its effect will, according 
to Saint-Venant’s principle, be large only near the ends 
of the beam, and will diminish rapidly with increasing 
distance from the ends of the beam. 

The second component stress distribution is now 
analyzed by the methods of the mathematical theory of 
elasticity for plane stress. 

The boundary conditions are as follows 


(1) até= +1, o,f, = 77,,= 0 
(2) at¢ oP 

oO = oO oO = Ou T = T () 
(3) at & O, 7 + | T = 2y71 


where op, and oo, are the negative of the longitudinal 
stresses in the flange and web, respectively, of the beam 
of infinite length. 

In an exact treatment of the problem, two additional 
conditions would have to be satisfied at the joint be- 
tween flange and web. The first of these is the con- 
tinuity condition that longitudinal strain in the flange 


and web should be equal at this joint. In terms ol 
stress this condition may be written as follows: 


(4) até =0,7 = +1, o. =<0,,.— vo, 


The second additional condition arises from the fact 
that the curve of displacement of an orginally straight 
line in the x-direction in the flange is discontinuous 10 
slope at x = 0, due to the shear flow into the flange from 
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stress measurements 


Experimental setup for I-beam heating and therma 
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the web. Because of this discontinuity, it is convenient 


to treat the two halves of the flange lying in the positive 


and negative directions of x separately, and to define 
a condition establishing continuity between them. 
Such a condition is one requiring that, along a line in 
the flange coincident with the joint, 


Ou; 3 = () 


where w, represents displacement in the x-direction. 
This condition may be expressed in terms of stresses 


as follows: 


(5) Até == (), (mu ¥) (O72, Or) = Oo. O§& 


- 


An exact determination of the stress distribution is 
not feasible in the present case. This distribution may, 
however, be determined approximately by application 
of the principle of least work in conjunction with the 
Ritz procedure of representing the stress function for 
web and flange regions as a superposition of a finite 
number of assumed functions in initially unknown pro- 
portions. '* 

It was found, after considerable investigation, that 
assumed functions in the form of polynomial expres- 
sions would be most suitable for the satisfaction of 
boundary conditions and convenient in subsequent 
mathematical manipulations. In consequence, the 
following series expressions for the stress functions in 
the flange and web regions were chosen: 


p p 
= b, + DY Anmnl— — 2G? — 1)72E"E"— (15) 
0 


m=On 


where 


gos = LS S o0,dé dé 


Dp 
Qu = Dow - > : ® 5,9" = 1) (¢* a, Ss ha aa (16) 


m=0 n=0 


PSS cogdn dy 


It should be noted that the expression for @¢,, as 
written, applies only in the range of positive £ This is 
unavoidable for reasons discussed previously in connec- 
tion with the development of boundary condition (5). 
Because of symmetry, the strain energy in the two 
halves of the flange will be equal, so that only one side 
need be considered in the energy minimization process. 

Continuity between the two halves could be estab- 
lished on an exact basis by the application of boundary 
It was found, however, that an exact 


where Dow 


condition (5). 
satisfaction of boundary conditions (4) and (5) greatly 
increases the difficulty of obtaining a satisfactory solu- 
tion with a reasonable number of assumed functions. 
Even if boundary conditions (4) and (5) are not ap- 
plied explicitly in the analysis, the energy minimi- 
zation process will tend toward their satisfaction. 
This follows from the fact that minimization of 
the strain energy on an exact basis leads to the 
true stress distribution only if the forces applied at 
the boundaries of the various regions do not change 
during the variation process or if no net work is done 
by changes in these forces (see reference 12, p. 166). 
The first of these conditions is seen to be automatically 
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satisfied at the free edges of the beam by each of the 
At the joint 
between flange and web the shear forces will be equal in 


assumed stress functions individually. 


magnitude and opposite in direction in each because oj 
boundary condition (3). If changes in these shear 
forces are to do no net work, the displacements through 
which they act must be equal. Consequently, this 


provides for the satisfaction of boundary condition 


(4). The assumed functions provide for stresses 
in the flange at its centerline (£ = O), and, conse 


quently, changes in these stresses will do no work only 
if the displacement through which they act is zero 
This, in turn, provides for the satisfaction of boundary 
condition (5). 

On the basis of these considerations, the analysis 
was carried out with boundary conditions (4) and (5 
eliminated from explicit consideration. Boundary 
conditions (1) and (2) are seen to be satisfied by the 
assumed stress functions. Satisfaction of boundary 
condition (3) leads to the following relations between the 
unknown coefficients A,,, and B,,,: 


n 


(y/w)(2A9, — Ay) = — » is Dp. tf =e 1,2. t 
” 0 


17 


Eq. (17) provides p + 1 relations which permit the 
expression of » + 1| of the unknowns as linear functions 
of the remaining unknowns in the total set of 2(p + 1) 
unknowns. 

The strain energy in the beam may be expressed as 


follows (see reference 12, p. 167): 


*! | 
V = (401 6,/E) | jor? + o,7 — 2voz,0,, 4 
/7V /7 0 


71 fl 
" a? — . ( 9 
2(1 + v)t22}dé d¢ + (2s/6,/E) | | + Oy_* + 


/7 0 0 


C<.* — Qe, s- 
w u u 


+ 2(1 + v)ry2,2}dn dé (18 
With 


or, = 0°,/0z", 0, = 0°,/Ox*,  7z,, = — (0°, 0x02 


and corresponding relations for the web region, Eqs. 
(15) and (16) may be introduced into Eq. (18) along 
with Eq. (17) to yield a quadratic form in the independ- 
ent unknown coefficients plus a sum of linear terms in 
these unknowns. 

The strain energy is then minimized by equating the 
partial derivative of the resulting expression with re- 
spect to each of the independent unknowns in turn to 
zero. ‘This yields a set of (p + 1)(2p + 1) simultane- 
ous linear algebraic equations in these unknowns 
which may be solved to yield the information needed 
in the determination of the stress distribution in the 
beam. 

The preceding analysis was applied to the case, ¥ = 5, 
w = 2, y = 1. The loading corresponded to a condt- 
tion of thermal stress resulting from transient heating 
as described in Section (4), with AX = 2.42, 7, = 0 
These parametric values were chosen to correspond 
with those appropriate to the experiment described in 
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SIMPLE 


the following section. The temperature distribution 


was that corresponding to the maximum temperature 
differential between the edges of the flange and the 
middle of the web. The stress distribution for the case 
of an infinitely long beam was first determined by the 
method of Section (5 This stress distribution was 
then reversed in sign and applied as an external loading 
on the ends of the finite beam to determine the end ef 


\ total of eighteen terms in the two stress func 


fect 
tion expressions were employed; that is, m = 0, 1, 2, 
n B12. 


[he distribution of normal stress in the middle cross 
section of the beam and the distribution of shear stress 
in the web-flange joint are shown in Figs. 14 and 15, 
respectively. The former distribution was found to be 
only slightly modified by the effect of the free ends. 
The waviness of the latter distribution in the central 
portion of the beam is probably not realistic, but is 
rather a consequence of the limitations of the approxi 
mate solution. It is believed, however, as a result of 
other computations performed, that the analysis pro 
vides a good indication of the magnitude of the maxi 


mum shear stress. 


7) EXPERIMENTAL INVESTIGATION OF THE TRANSIENT 
TEMPERATURE DISTRIBUTION AND THERMAL STRESSES 
IN A WIDE-FLANGED I-BEAM 


The theory of the preceding sections was checked 
experimentally by means of the apparatus shown in Fig. 
; inch steel 


12. The specimen used was built up of 
plate, with the web-flange joints brazed in order to 
provide an unimpeded path for heat flow from the 
flanges into the web. The remaining dimensions of the 
beam were 2/ = 30 in., 2b = 10in., 2s = 5 in. 

Convective heat transfer into the outer faces of the 
flanges was simulated by the transfer of heat from large 
heated aluminum blocks across a small air gap (0.007 
in.). The proper gap was maintained by means of 
suitably located spacers, and the effect of flange war- 
page was eliminated by providing an initial concave 
outward curvature and exerting enough pressure across 
the two heating blocks to bend the flanges back to a 
flat shape. 

The blocks were preheated by means of electric re- 
sistance rods, and then brought up to the specimen and 
clamped tightly. The surface temperature of the 
blocks was intended to simulate an environment tem- 
perature varying as a step function, and the thermal 
conductance of the air gap was intended to simulate 
the convective heat transfer coefficient. As heat 
flowed across the gap, there was found to be a drop in 
the surface temperature of the block, thus violating 
the assumption of a step function variation of environ- 
ment temperature. However, there was simultane- 
ously an increase in the conductivity of the air in the 
gap due to a rise in its temperature, and the two effects 
were found to be largely compensating. The value of 
the parameter \ was estimated to be 2.42. 
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Temperature was measured at selected stations on 
the inner faces of the flanges and on the web surfaces 
by means of electric resistance gages and stress was 
measured at similar locations by temperature-compen 
sating wire resistance strain gages. Both temperature 
and stress were recorded continuously by means of os 
cillograph equipment. 

The temperature and longitudinal stress distributions 
in the middle cross section at the time of maximum 
temperature differential are shown in Figs. 13 and 14 
respectively. 

Shear stress was measured by means of strain rosettes 
at locations on the surface of the web near one end of the 
beam and close to the web-flange joint. The experi 
mental values obtained at the time of maximum tem 
perature differential are shown in Fig. 15 along with 
the corresponding theoretical distribution as obtained 
from the analysis of Section (6 The experimental 
results indicate that the peak value of shear stress 
occurs closer to the end of the beam than is indicated 
by the approximate theoretical result, and is slightly 
lower in value. 

The results in general indicate good agreement be 
tween theory and experiment, and encourage confidence 
in the theoretical method. 
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The Effect of Amplitude on the Torsional 
Vibrations of Solid Wings Subjected to 
Aerodynamic Heating’ 


JOSEF SINGER* 
Polytechnic Institute of Brooklyn 


SUMMARY 


The torsional vibrations of a long thin solid wing subjected to 
aerodynamic heating are analyzed by the Rayleigh-Ritz method. 
The effect of the change in thermal stresses due to finite ampli- 
tudes is included. To simplify the solution of the resulting non- 
linear equations, a single-term deflection function is assumed, and 
the assumption is justified. An explicit expression is obtained for 
the relative frequency for any given temperature distribution and 
initial amplitude of vibration. A typical example is worked out 
to indicate the significance of the alleviating effect of the initial 
amplitude. 

The assumption of harmonic vibrations used in the Rayleigh 
Ritz approach is examined by a more accurate analysis yielding 
elliptic time functions and is shown to be satisfactory for practical 


purposes. 


SYMBOLS 


Fixe coefficient of amplitude function 

Ayn» = coefficient of main term 

a = semispan of wing, in. 

b = semichord of wing, in. 

ic. = effective polar moment of the cross section, in.4 

D = bending stiffness = FA3/12(1 — v?), lb. in 

E = Young’s modulus, psi 

F = time function 

G = shear modulus, psi 

g = gravitational constant, in. per sec.? 

h = thickness of wing, in. 

H = second moment of area about vertical axis of symmetry, 
in. 

I second moment of the fictitious edge force ha,/R, in.4 

K = constant defined by Eq. (9), in.® 

k = twist per unit length, rad. per in 

m = number of nodal lines parallel to x axis 

n = number of nodal lines parallel to y axis 

R = coefficient of thermal stress = aET), psi 

7. = temperature, °F. 

T) = temperature parameter 

V = potential energy, lb. in. 

w = deflection, in. 

a = coefficient of thermal expansion, in. per in. F 

+ = specific weight, lbs. per in.* 

r = parameter defined by Eq. (17) 

Mm = parameter defined by Eq. (18) 
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v = Poisson’s ratio 
' 

w = circular frequency of vibration, rad. per sec 

wo = circular frequency of vibration when 7; = 0, rad. 


sec. 
DISCUSSION 


HE RAYLEIGH-RITZ METHOD is often used to obtai 
ye approximate solution to the problem of the vibra 
tion of a rectangular plate, all edges of which are fre 
In this method, the maximum potential energy ; 
equated to the maximum kinetic energy for harmon 
vibrations, yielding an expression for the angular fr 
quency in terms of the amplitude w(x, vy). The ampli 
tude function w is taken as a series 


Zz > A pq? p(X) Gq(y) 
q=0 


p 0 @= 


and the coefficients 1,, are adjusted to minimize th 
expression for the angular frequency. 
The method can also be used for torsional vibrations 





of a thin solid wing subjected to aerodynamic heating 
The resulting thermal stresses reduce the torsiona 
rigidity and, hence, the frequency of vibration. 
However, since the thermal stresses change with ampli 
tude,’ the amplitude of vibration affects the fr 
quency—i.e., the larger the amplitude the less the 
duction in frequency due to the thermal stresses. Thi 
alleviating effect is noticeable even for small vibrations 

For simplicity, a long wing is considered for which th 
local stress effects at the tips can be disregarded and tht 
thermal stresses can be approximated by those of a wing | 
of infinite aspect ratio. Hence, as in reference 1, tht 
thermal stresses are o, = R f (y) and o, = 7,, = 0 whe 
the wing is in the neutral position (see Fig. 1). | 

As the wing twists, displacements w perpendicular t | 
the plane of the wing give rise to a strain 
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Fic. 1. Solid wing and thermal stress distribution in its mid 
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(1/2)(Ow/Ox)* — 


/ b eb 
(: / h ay) | (h/2)(Ow/Ox)*dy (1) 
« b e b 


for a long wing with no force or moment resultants ap- 


plied to its ends. The effective thermal stress is, there- 


fore, 


o, = Rf (y) — (E/2) X 


/ feb ah 
c Ox)? — (; | h dy) | h(Ow /Ox)*dy (2) 
7 b e b 


compressive stress being considered positive.’ 
The maximum potential energy becomes 


(1/2) | | | D(V2w)? — 2(1 — v)[(O°%w/Ox?) X 
e be r 


| | h xX 
e be a 


(Ow/Ox)? — 


(0°w/Oy?) — (O°w/Oxdy)?|idxdy — (1/2) 


; R f (y)(Ow/Ox)? — (1/2) (7/2) 


| gb eb 
( | h dy) | h(ow Ox) %dy | Ou dx)? baxdy (3) 
ny b e b 


Note that the potential energy of the change in thermal 
stress, the second term in the second integral of Eq. (3), 
is multiplied by an additional factor (1/2). This factor 
appears since the change in thermal stress varies from 
zero to its value at the amplitude w while Xf (y) remains 
constant. 

The circular frequency is then 


V/ 1 (y/g) | | h widxdy | 
e be a 


If the assumed amplitude function has more than one 
term, a set of nonlinear simultaneous equations are ob- 
The solution of these equa- 


(4) 


~ = YD 
Ww — 


tained after minimization. 
tions may present considerable difficulties. 
Pavlik® while investigating the 


However, found, 
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vibrations of rectangular plates of constant thickness 
with all edges free by the Rayleigh-Ritz method, that 
the coefficients of all terms of the amplitude function 
other than the main term are small and can be neglected 
for approximate frequency calculations. The main term 
iS <1 mnIm(X)¢n(y), Where m and nm are the number of 
nodal lines parallel to the x and y axes for the mode of 
vibration considered. 

This, and other similar results (some also verified ex- 
perimentally —e.g., reference 7) led to Warburton’s ap- 
proximate frequency formula for rectangular plates of 
constant thickness with various types of boundary con- 
ditions,* which is based on a single-term amplitude func- 
tion and has been shown to be in good agreement with 
more exact calculations. 

For the first torsional mode, m = | and n = 1, the 
single-term deflection function (or main term of the 


series) is 
w= kxy (5) 


where & is the angle of twist per unit length. 

Pavlik and Warburton have justified this form of de- 
flection function only for rectangular plates of constant 
thickness without any axial loads. For the case of self 
equilibrating axial stresses of the type resulting from 
aerodynamic heating the author has compared the fre- 
quencies calculated for the single-term deflection func 
tion of Eq. (5) with those calculated for a three-term 


deflection function 


w= Ayxy + Agx*y + Aisry® (6) 


He found errors of the order of 5 per cent near the 
critical temperatures where the axial stresses are larg 
est. Errors of the same magnitude are found for plates 
of parabolic chordwise thickness variation. 

These results of small deflection analysis justify the 
assumption that, for the purpose of investigating the 
effect of amplitude on the vibrations, Eq. (5) represents 
a first approximation also for a thin wing of chordwise 
varying thickness subjected to thermal stresses given by 
Eq. (2). 

With this assumption, Eq. (4) becomes 


w? = [GC — RI + (E/2)K k?|/[ya*H/3g| (7) 
where GC is the torsional rigidity and 
b 
[ = [ h y*f(y)dy (8) 
J» 
eh / rb eb 7 
K = h E —y (; / h dy) | hy? dy dy (9) 
e b le b e t J 
+b 
H = hy*dy (10) 
« b 


Hence, the frequency depends on k, the twist per unit 
length. Since the amplitude of vibration depends on the 
initial amplitude, the frequency can be estimated 
for any initial k from Eq. (7) if damping is neglected. 
The relative frequency, the ratio of the frequency to 


that for no temperature variation, is given by 
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(w/w)? = 1 — [RI — (E/2)K k?]/GC (11) 


Fig. 2 presents numerical values of this ratio com- 


puted for a wing of constant thickness h = 1/4 in. (flat 
plate) when the chord 26 is 20 in. The temperature dis- 
tribution is given by 7 = 7) + 7\(y/b)*; the material 


constants are E = 29 X 10 psi, a = 7.5 K 10-®°F.—!, 
and Poisson’s ratio v = 0.5. The effect of the initial 
amplitude on the relative frequency is significant even 
for small amplitudes—e.g., for k = 0.0020 (correspond- 
ing to a maximum tip rotation of 4.1° for a wing of span 
SO in.), the relative frequency is 0.471 at the critical 
temperature instead of zero as given by the small de- 
flection analysis. 

In the Rayleigh-Ritz approach used above the vibra- 
tions were assumed to be harmonic. Actually, the time 
functions become elliptic functions when the changes in 
thermal stresses are taken into account. This can be 
shown by assuming, as before, that Eq. (5) represents 
the deflection for the first torsional mode of vibration. 

‘Hence, for a wing constrained by the above assump- 
tion to vibrate in this mode, the displacement at any 
time is 


w= kxyF (t) (12) 
where F is the time function associated with the mode 
kxy. Substitution of Eq. (12) into Eq. (3) gives the 
potential energy* 


V = GCak?F? — RI ak?F? + (E/4)aR‘KF! (13) 


and the kinetic energy 


b a 
T = (¥ 2) f f hwdxdy = ya*k*HF?/3g (14) 
—b a 
where /, K, and // are given by Eqs. (8)—(10). 
Lagrange’s equation becomes 
(d/dt)(OT/OF) + dV /OF = 0 (15) 


since F can be considered the generalized coordinate. 
If the operations indicated are carried out, one ob- 


tains 
F+)\F+ uF? = 0 (16) 

where 
A = (GC — RI)/(ya?H/3g) (17) 
w= [(£/2)K k?]/(ya?H/3g) (18) 


* V in Eq. (8) is the maximum potential energy only when w is 
the maximum displacement—i.e., the amplitude—as in the Ray- 


leigh method. 
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Eq. (16) is the same nonlinear ordinary differen 
equation as the one obtained by Burgreen’ for the f; 
vibrations of a pin-ended column with constant distay 


between pin ends. The initial conditions 





‘ 


w(x,v,0) = 0 


w(x, y,0) = kxy and 
and the boundary conditions 


F(O) = 1 and F(0) =0 


are also the same. Burgreen’s solutions for F in terms 


Jacobi cosine and delta cosine elliptic functions : —————— 


therefore for the present problems also, the form of th 
solution depending on the value of Xu. 

However, it can be shown that, except for near} 
vanishing frequencies, the shape of the time function 
fairly close to that of the circular cosine function a 
that the relative frequencies obtained by the Rayleig} 
Ritz method differ only slightly from those obtained | 
the more accurate analysis unless k becomes large 

For the example considered, this can be seen fro: 
Fig. 2 where the more accurate relative frequencies a 
also plotted for comparison. In view of the approx 
mate nature of the physical assumptions, the Rayleig! 
Ritz method appears to be sufficiently accurate for th 


present purpose. 
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Laminar Boundary Layer With Nonadiabatic 
Rotational Free Stream 


Franklin D. Hains 

Aerodynamicist, Research Division, Bell Aircraft Corporation, 
Buffalo 

April 5, 1957 


ob POSSIBLE IMPORTANCE of vorticity behind curved shocks 
in high-speed flow was discussed by Ferri.'! Recently, shear 
flow over a flat plate was considered by several authors.? 

This note describes a method of solving the flow over a flat 
plate with heat conduction and transfer of momentum into the 
boundary layer from the external flow. It is assumed that the 
normal gradients of temperature and velocity in the free stream 
are small so that a boundary layer exists adjacent to the plate 
where these gradients are large. 

The boundary-layer equations in the Crocco form are 


2 


T°T,, — ppur, = 0 (1) 


(1 — P)rrh, + (h,, + P)r’? puuPh, = 0 (2) 


Prandtl Number and linear viscosity-temperature 


assumed 


A constant 
relation are 


The boundary conditions are 
h=h, t 2 
T = Uo, h, = B/a { 


u 


where a and £ are the slopes of the velocity and enthalpy profiles 
at the outer edge of the boundary layer 
A solution to Eqs. (1) and (2) is assumed of the form 


Z fi(u) b*-"/+/2x $) 





i 
i 1 
. k x 
h= ye t(u)vi-! 4 } 2 a,x"¢,,(u)v'~! 5) 
1 1 n O1 l 
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where 


® = (a’v,,x/u,,%)! 2) 
{ (6) 


VY = Bu? ah, 


It is assumed that © and W are small so that higher order terms in 
# can be neglected. The constants a, in Eq. (5) are determined 
by the wall-temperature distribution which is represented by a 
polynomial in +/ x. 

Substitution of Eqs. (4) and (5) into Eqs. (1) and (2) and equat 
ing to zero coefficients of like powers of ® yields 


fifi "4 ppu = 0 (7) 
fife” + 2fi"fe = O (8) 
(1 — P)fi't’ +fili” + P) =0 (9) 
(1 — P)[fifi'te’0 + (fife) ’tr’] + 2fifetr” 4 
f7t2"0 + 2fifoP — puuPét, = 0 (10) 
(1 P\fifi'gin’ + frei,” — 2ppuPng,, = 0 (11) 
(1 — P)[(fif2) fin’ + Sift 'gon’O) + fr?gen”@ + 
ZAifoein” — (2n + 1) ppudPg,, = 0 (12) 
where 0 = Bu,,/ah, (13) 
The boundary conditions are: 
u=0: fi’ =f.’ =0 ) 
u=Uu, fi =t = bt’ = gr, £2, Lin. go,’ = 14) 
h=h,, t’ = B/a, fe = Vo Patten? 


The unknown quantities in Eqs. (4) and (5) are obtained from 
the solutions of Eqs. (7) through (12). Eq. (7) can be transformed 
into the Blasius equation, while the remaining equations are 
linear. Irregular singular points occur in the remaining equations 
at u = Ue so that asymptotic solutions are assumed near the 
outer edge of the boundary layer. 

At the present time, numerical solutions are being obtained 
for problems of aerodynamic interest 
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Extension of ‘‘An Aeroelastic Parameter for 
Estimation of the Effects of Flexibility on the 
Lateral Stability and Control of Aircraft’’! 


William P. Rodden* 
Senior Flutter Engineer, North American Aviation, Inc., Los Angeles 
April 1, 1957 


A SIMILARITY PARAMETER is suggested in reference 1 to simplify 
the modification of rigid stability derivatives in a flexible 
stability and control analysis. It is the purpose of this note to 
make a refinement of one of the approximations and to present 


some experimental correlation of the method 


* The author wishes to thank the Douglas Aircraft Company, Inc., in 
particular Kermit E. Van Every and Ervin R. Heald of the El Segundo 
Division, for permission to use the data of Fig. |] 
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Fic. 1. Experimental correlation. 








We write, as before, the twist at the wing tip, e, as an infi 
series of twists: 
€ = € + €& + 6 4 
where e, is the twist resulting from the initially applied loads 
€1 = €i)a + €p(pb/2V) + ep 
and each successive €;,; is caused by the loads resulting from, 
Now we expand each ¢ in terms of its components from t 


initial value 


E(i41)5/ Cig = ks 
E(i+1)p/ Fip = k, 
€ 1 €i, = k; 


The series for « becomes a geometric series whose sum is 


e = [(€5652.)/(1 + &s)] 4 [ep( pb 2V)/(1 + &,)] 4 

€-p (1 +4 k 5) 
provided that k > I As pointed out recently by Gaugh a 
Slap,? the condition k < —1 corresponds to antisymmetrical win 


divergence. They also have shown in their earlier work? that E 
(6) is still valid for k > 1 by extending the radius of converge 
of the geometric series. The modified rolling-moment equati 


now becomes 


C'1,p = C'1,5.4 + C',,(pb/2V 
where 
Cp, = Cy 11/(1 + k;) eS 8 
Cy = Ci + [1/0 ks) (Cr) 565 
C',, = C, + [1/01 kp)| (Ci) re l 


in the notation of reference 1, and the various k become 


bk, = —(de/dL4) ;(Cy,) ;qSb 
ks = (0e/OLe)5 C1.) aqSb 
k, = (6¢ OL) p( Ci) pqSb l 


where the twist aeroelastic parameters, Je/OL,, now are fou 
from each twisting mode. This extension from the developmet! 
given previously eliminates the restriction of computing & from 
mean twisting mode and increases the generality of the method 
It has been suggested? that a better example for correlatio 
would be one with a lower stiffness level than that previously pr 
sented. In Fig. 1 an experimental correlation of aileron effectiv' 
ness is shown based on some low-speed wind-tunnel tests of 
flexible swept-wing model by the El Segundo Division of Douglas 
Aircraft Company, Inc. The estimate was based on the ailero! 
effectiveness obtained from rigid models and the structural 
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| A Short Static-Pressure Probe Independent 


mads 


from ‘ 


‘om t 


Number 








pressure error in per cent of impact pressure as in reference 2 


Erik L. Mollo-Christensen, * 
Martuccelli 
Massachusetts Institute of Technology, Cambridge, Mass 


April 8, 1957 


of the boundary layer are insignificant 


I REFERENCE 2, Korkegi presents means of static-pressure 
error compensation by suitably shaping the boundary ii 


are thought to be of general interest. 
terms of per cent static-pressure error rather than in terms of 


of Mach Number? 


Marten T. Landahl,** and John R 


*** 


SUMMARY 


Using slender-body theory it is possible to shape the nose of a pressure 
probe such that it measures stream static pressure independent of Mach 
The pressure holes can be placed so far forward that the effects 


Such a probe has been built and 


tested, and the probe performs as predicted by theory 


INTRODUCTION 


the 


vicinity of the pressure orifice 
rhis note provides some additional experimental results which 


The results given are in 


») 


DISCUSSION 


Consider the problem of shaping the nose of a body of revolu 
. tion such that, at a distance x from the nose, the pressure equals 


static pressure in the stream far ahead of the body for any free 


stream Mach Number 


Slender-body theory gives the following formula for the pres 


sure coefficient :! 


C(x) = (p(x) 


| where 


p 
p 
I: (1/2)pl? 
\ 
R(x 
S(x 
V/ 


** Research Engineer 








S’(0) In x 4 


Pol/(1/2)pU? = 
S"(x) In [2/R(x)V M? 1] + 


. 


. / 
| In (x t) dS"(t)¢ [R’%(x)}? (1) 
e tT 


local pressure 


= pressure in the undisturbed stream 


= dynamic pressure 


Inspecting Eq. (1), it can be 
pressure holes will be drilled in the probe, S”"(x1) = 0, Cp(x1) be 


comes independent of ./ 


+ This work was in part supported by the U.S 
nautics, under Contract No 


* Assistant Professor 


*** Research Engineer 


distance from the nose of the body 

body radius at x 

m|R(x)}2, body cross-sectional area 

Mach Number of the undisturbed stream 


seen that if, at x,, where the 


Therefore, impose the condition 
S"(x,) = 0 (2) 
Navy, Bureau of Aero 
NOas 56-556-d 
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The other terms in C,(x,) can be made to vanish by a suit 
able choice of S(xv) over the interval 0 < x < x 

For example, S(x) can be expressed as a fourth-degree poly 
nomial in x and the coefficients chosen so as to make C,(x)) 0 


The result is 


R(x) = (6/V w)V x? 0.968x* + O.317x4 (3) 
where 6 is a thickness parameter which can be chosen freely 

The resulting probe contour is shown in Fig. 1 

If now three holes are drilled at x,, 120° apart around the cir 


theory shows that the 


cumference of the probe, slender-body 
mean pressure in the three holes when the probe is at 
of yaw is the same as when the probe is aligned with the stream 


an angle 


A probe of the dimensions shown in Fig. 1 was built and tested 
in the M.I.T. 2 X 2'/2-in. supersonic tunnel. The probe pressure 
was compared with the pressures measured in static orifices in a 
flat plate in the wind-tunnel wall. The wall boundary layer was 
removed at the leading edge of the plate. Visual inspection of the 
flow using schlieren apparatus provided information of probe 
wall interference, and such interference does not occur in the 
results shown in Fig. 2, where | p(x) P.|/P. is plotted versus 
Mach Number 

The results show an almost the Mach 
This error is probably due to imperfections 


constant error over 
Number range 
in the probe contour which is difficult to machine accurately for 
the size used. With some further refinement, such a probe 


may prove a useful tool in wind-tunnel work and in flig ht testing 
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On the Exact Solution of Compressible Couette 
Flow With Mass Addition and Binary 
Diffusion 


Nathan Ness* 
Research Scientist, Missile and Ordnance Systems Department, 
General Electric Company, Philadelphia, Pa. 


April 8, 1957 


SUMMARY 


An exact, closed-form solution of the system of equations describing 
compressible Couette flow with mass injection and binary diffusion is ob 
tained. Velocity, enthalpy, and diffusion profiles are expressed in terms 
of simple integrals for the completely general variation of all fluid properties 


SYMBOLS 


a ue?/(e "40 — 1)2 

( mass fraction of species i; C; pi/p,# = 1,2 
specific heat at constant pressure of mixture 
specific heat at constant pressure of species 7; 7 ‘2 
Die binary mixture diffusion coefficient 

I definite integral defined by Eq. (23) 


/ static enthalpy, hi Cp aT 


K constant of integration, defined by Eq. (21) 
coefficient of thermal conductivity of mixture 
Fy Prandtl Number, Pr; MC p,/k 
S Schmidt Number, Sc u/pDie 
1 absolute temperature of mixture 
t velocity of mixture parallel to wall 
velocity of mixture normal to wall 
absolute velocity of species 2 normal to wall 


y coordinate normal to wall 
” ey 
8 Pridn pv Pri(dy/p) 
0 0 

\ film thickness 

. y 

pv (dy ‘) 

0 
yu coefficient of viscosity of mixture 
p density of mixture 


ry 
¢ pv Sc(dy/p) 
0 


Subscripts 
l injected fluid 
2 mainstream fluid 
D diffusion film layer 
outer edge of film layer 
] temperature film layer 
velocity film layer 


at the stationary wall 


DISCUSSION 


LLINGWORTH! obtained an exact solution for the system of 

equations describing compressible Couette flow; Liepmann 
and Bleviss? extended the analysis to include the effects of disso- 
The 


with 


ciation and ionization. present analysis considers com 


pressible Couette flow mass addition and diffusion. An 


exact solution results with completely general variation of all 
fluid properties. 
mixture, com 


The system of equations describing a binary 


pressible Couette flow is® 


Continuity : (d/dy) (pv) = 0 l 


Momentum: pu(du/dy) (d/dy) |u(du/dy)} 2 
The author would like to express his gratitude to Drs. Eckert, Hartnett 
and coworkers at the University of Minnesota for helpful discussions during 


the prepacation of this note 
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Diffusion : pu(dC,/dy) = (d/dy) [pD2(dC,/dy)} 


Energy: 


puc,(dT/dy) = (d/dy) [k(dT/dy)| u(du/dy)? 4+ 





pio Cp, Cn) LET dy) (dC, /dy 


The boundary conditions are: 


y= 0: a=Q0, pv =(pv).,, T=Te, GC c 


y = A: umu, T=, G=C, 


The solution of Eq. (1) is py = constant. The boundary 


dition y = 0, pu = (pv)» indicates that the flux of the mix 
normal to the main flow is constant throughout the film | 
Consequently, to satisfy continuity across the flow, it is necessy 
that the moving wall be considered porous with a suctio 
equal to (pv) ,». 

Defining a new variable n by 


ng v 
n = pi | (dy/p) 
* 0 ‘ 
enables the momentum equation to be written in the form 


du/dn = (d/dn) (du/dn) 
whose solution, satisfying the boundary conditions y = 0 
0O,y = A,,u = 4u, is 
1)/(e7Ar 1) 
In a similar manner, the introduction of a new variable ¢ 


fined by 


g = pl | Sc(dy/p) 
0 
enables the diffusion equation to be expressed as 
dC\/dg = (d/d¢g) (dC,\/de) 


whose solution, satisfying the boundary conditions y = 0, ( 
Cis, 7 = Ap, C = ti is 


C;) (Ce Ci, = (e* 1)/(e¥ Ap 


) 


( Cu 


The absolute velocity of the air is equal to its diffusion v 


locity plus the observed macroscopic velocity of the mixture 
The diffusion velocity is given by Fick’s law; hence, for tl 
y-direction, the above statement gives 
= (Dio/C2)+(dC/dy) + 12 
or, since C, = p»/p, the resultant flux of the air normal t: 
walls is 
py» p Dy d Ce dy) + ps 


The term pDyw(dC2/dy) represents the mass flow of air by dif 
fusion across any control surface parallel to the walls; the terr 


pw represents the convective mass flow of air across the 


control surface. To ensure continuity it is necessary at 
point in the flow that 
pw = pD(dC2/dy) 
Therefore, from a comparison of Eqs. (13) and (14), it follow 


that pve = O everywhere in the flow The implication of tl 
result may be more readily seen by considering the relationshy 
that exists between the flux of the mixture and the flux of tl 


respective components: 
pe pil T pws 


As pv. = 0, it follows that the flux of the mixture across t 
flow is equal to the flux of the coolant in the same directior 
Using the condition that the sum of the mass fractions of th 
air is unity—i.e., C, + C, = 1—enables Ee 


coolant and the 


(14) to be written 


pv = [pDye/(1 C,)] -(dC,/dy) 


Eq. (15) is substituted into Eq. (4) to give 


pvcp(dT/dy) = (d/dy) |k(d7T /dy u(du/dy)? 


pr(1 ( (Cp, Cp) (dT /dy 


Since Cp = 


pvc p,\ 


Defining the 


nables the la 
d 


Using Eq. (9) 


Introducin 


gives Eq. ( 19 


dh, dp 


Integratin; 


(dh 
where Ky is 


ondition 


The const: 


e enthalp 
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Cyep, + Coep., the energy equation simplifies to 


Since Cp 


pcp (dT /dy) = (d/dy) |k(dT/dy)] + w(du/dy)? (16) 


nefining the Prendt! Number by Pr; = wcp,/k and the enthalpy 


0 


e7 
hy -{ cp dT, dhy = cp, a7 (17 


nables the last equation to be written in the form 


dh,/dy = (d/dn) (dhi/Pridn) + (du/dn)? (18) 
Using Eq. (9), Eq. (18) becomes 
dh,/dyn = (d/dn) (dhi/Pridn) + a e*” (19 
vhere a = u,*/(e"A 1)? (20 
Introducing a new variable 6 defined by 
ey 
B = pv | Pri(dy bw) (21) 
* 0 


gives Eq. (19) in the form 


dh,/dB = (d/dB) (dh,/d8) + (a 


Integrating once results in 


(dh, /dB) h, = (a/2e K 


where K, is the constant of integration. This is a first-order, 
linear differential equation whose solution satisfying the boundary 


ondition 


(a/?2) | e 9 dp + Kye p 1) 


The constant A, is evaluated from the boundary condition 


v Ar, 1 ls “ hy = hn, 


dp 


‘ 7 3 j 
2(h hy, + [(e? | ePA L)Je%Aa 
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A Note on the Sonic-Wedge Leading-Edge 
Approximation in Hypersonic Flow 


M. H. Bertram and D. L. Baradel! 
Langley Aeronautical Laboratory 
April 18, 1957 


NACA, Langley Field, Va 


oe in the Princeton helium tunnel! have cast doubt 


on the validity of approximating the blunt leading edge on 


flat plate by a sonic wedge as proposed by the senior author 


Recent theoretical papers*: * using a blast-wave analogy have 


indicated the parameters governing the effect of a blunt leading 


edge on the pressures and shock shape over a flat plate, at least 


in the vicinity of the leading edg« The parameter governing 
the bluntness effect on surface pressure is A Vo *kt/x, where 
k is the nose-drag coefficient and ¢ is the nose thickness. The 
corresponding shock shape is a function only of V/k(x/t Phe 


present authors have re-examined the sonic-wedge leading-edge 


concept in the light of these parameters 
la are the surface pressures given by rotational 
7/5. Also in 


a cylindrical leading-edge flat plate 


Shown in Fig 
characteristics theory,® this Figure is the 
characteristics solution for 
By Newtonian 


hypersonic approximation, the equation given by 


theory for a cylinder, k (2/3)C; To 


Lees® may be 
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2 r 
IO Le eS Saye ee: 


used: 


Cp mez = [ly + 3)/(y + 1] {1 — 2/f(7 + 3)Mo%4} 0 


“as 


k= (2/4) {1 — (y + 3)/[(y + 1I)Ma3)} (2 


The correlation of the various characteristics solutions in Fig 
la is good except very near the shoulder, and there is reasonable 
agreement of the correlated portions of the characteristic solutions 
with the blast-wave analogy toa value of A of about 100. There 
is a compression above the pressure obtained immediately after 
the expansion around the shoulder at each Mach Number shown 
\ quick first approximation to the pressure distribution for the 
sonic-wedge leading-edge flat plate (and probably the evlindrical 
leading-edge plate as well) is obtained by taking the pressure as 
constant at the shoulder value to the intersection with the cor 
relation curve or the line for the blast-wave analogy, which is 
used thereafter in its range of validity. A good fit to the cor 
relation curve is obtained by adding an empirical constant to the 
first-order blast-wave theory result, thus obtaining an equation 
analogous in form to Sakurai’s second-order three-dimensional 


result.’ In this case using the Cheng-Pallone equation,’ 7 = 
7/5, Ap/po = 0.112 K2/3 — (0.2, where 0.20 may also be used 
for y= 5/3. 


For the sonic-wedge leading-edge flat plate the shoulder pres 
sure is, to hypersonic approximation, 


Ap/p2 = (yMo?/2)a[1 — (b/M«?) 


where, for 7 = 7/5, a = 0.10082 and b = 6.5, while, for 7 = 
5/3, a = 0.092819 and b = 5.85 

Shock shapes corresponding to these pressure distributions are 
shown in Fig. lb. Characteristics theory indicates the shock 
shape to be independent of 7 only in the immediate vicinity of 
the leading edge to \/x as low as 5 but with an increasing ex- 
tent of correlation as ./. increases 

The question arises as to how many approximations to each 
characteristic point are required for an arbitrary accuracy. For 
Mo = 6.86, the characteristics solution was carried out to one, 
As for all these 


two, and three approximations for 7 = 7/5. 
solutions, the expansion at the shoulder was carried out in 1/10 


steps for the first degree and 1° steps thereafter. The exter; 


x/t over which the solution is valid is clearly shown in Fig 


In this case to get good accuracy for the pressure distribution j 
x/t > 100, at least a three-approximation solution is requir 
There is little difference between the one-, two-, and thre 
approximation solutions for the corresponding shock shapes 
Fig. la the solutions were all cut off at what was considered th 
limit of good accuracy, and the solution for Ma = 6.86 was, 
tended by the three-approximation result. 

Given in Fig. 3 is a comparison with the Princeton heliy 
tunnel data.!- Only data considered to have a small to negligih 
viscous effect are included, and k = Cpmaz [Eq. (2)] is assume 
Considering the Mach Number gradient in the Princeton heljy 
nozzle, the agreement between characteristic theory (two 4 
proximations) and experiment in Fig. 3a is considered B00 
This characteristic solution does not agree with the one present 
in Fig. 7 of reference 1. The differences between these two so) 
tions are large and are apparently attributable to an oversj 
plification of the flow model in reference 1. Of interest js ; 
fact that the characteristic solution presented in Fig. 3a con 
lates (excepting the region near the shoulder) with the air resy} 
in Fig. la if corrected to y = 7/5 according to the constay 
given in reference 4. 

The experimental points for shock shape, shown in Fig. 3b, 
between blast-wave and characteristic theory. However, t 


effect of the nozzle pressure gradient is unknown, and in t 


region very near the nose the differing nose shapes and shoei 


origins must be important as regards shock shape 
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exten REFERENCES compares very accurately with that predicted by other surface 
in Fig Hammitt, A. G., and Bogdonoff, S. M., Hypersonic Studies of the Lead theories —e.g., Krienes’ solution; electromagnetic-analogy meth 
ution j » Bdge Effect on the Flow over a Flat Plate, Jet Propulsion, Vol. 26, No. 4 ods; three-quarter chord line solutions; formulas used in the 
requir yp. 241-246, April, —- ; past such as the R. T. Jones edge-velocity method for elliptic 

t “ Mitchel H., Viscous and Leading-Edge Thickness Effects on the . . . . , 
* Bertram, abit: luading, whic ™~ 6, lea = 2 - 2.5) 
id thre pressures on the Surface of a Flat Plate in Hypersonic Flow, Readers’ Forum vading, < hich for A }, leads to Cra 2rA (A P + 2.5); and 
ipes journal of the Aeronautical Sciences, Vol. 21, No. 6, pp. 430, 431, June an experimental formula for rectangular wings, which for 4 > 
red theef 1954 6, leads to Crag = 20A/(A + 3). For other than elliptic loading 
Lees, Lester, and Kubota, Toshi, /nviscid Hypersonic Flow Over Blunt 4 : ‘ a : 
rac « 4 as § — ) »)) 
ita \ysed Slender Bodies, Journal of the Aeronautical Sciences, Vol. 24, No. 3 Cra = 24A/j{A + 2(A + 4)/(A 4 2) (1+ 47 
» 195-202, March, 1957 cian . _ : 7 ; . 
1ere is given by line theory, gives acceptable enginee 4 
1 heliy “Cheng, H. K., and Pallone, A. J., Jnviscid Leading-Edge Effect in Hype Sis g ~ ; ray bs a — ap ae Eg 
egligit nic Flow, Readers’ Forum, Journal of the Aeronautical Sciences, Vol. 23 accuracy For the A = 6 rectangular wing, r > 0 2 Putting 
a ~ E xo.7, pp. 700-702, July, 1956. A = 6 inside the brackets gives directly the previous two formulas 
sum ) . : : : = as 
ss Ferri, Antonio, Application of the Method of Characteristics to Super Applying simple-sweep theory to the section coefficient leads 
1 heliy: nic Rotational Flow, NACA Rep. 841, 1956. to acceptable values of C for swent wi | . . tel 
« - ‘ “ rt = « ( swe Z Ss av < 5 “ = 
two ay Lees, Lester, ‘‘Hypersonic Flow,” Fifth International Aeronautical Con-~ = ae I . > Len é bi =e sitetiiie-: a ay 
d g00 ence, (Proceedings), IAS-RAeS, Los Angeles, June 20-23, 1955; Institute elliptic loading rhen substitute (8( La/k COS Ag) for Cy, and 
inal tthe Aeronautical Sciences, New York, 1956 (8A /x cos Ag) for A, where tan Ag = tan A;2/8 and where A,» is 
sente : : : 
ical the sweep angle of the mid-chord line of the wing 
WO Sol 
VErS} 
t ist + 
1 Corr 
Tesult —_ 
nstan} Ryle of Thumb Equation for Predicting 
Lifting-Surface-Theory Values of Lift 
3b, 
as a — Vortex-Ring Cascade and Turbulence 
, n #7 
a ee ; ; . Transition 
“i ITA, Centro Tecnico de Aeronautica, Sao José dos Campos, 
Shock Estado Sao Paulo, Brazil y OR 
Vv ri 
April 18, 1957 ivian en 
Applied Physics Laboratory, The Johns Hopkins University 
Nii ; ; Silver Spring, Md 
— fF" WINGS a function is derived that is almost as easy to re April 29, 1957 
member as the familiar equation of lifting-line theory 
_~f Compare (elliptic loading; Cra, lift-curve slope; A, aspect rati . 
HT sac ! B Le I J Pore r | “MHE FACT that vortex rings tend to decay by a cascade process 
. Cha, = 297A/(A + 2) has been known for some time (see, for example, Sadron'), 
o = 27A/{A + 2[(A + 4)/(A + DI} and interest in the process has recently been revived by Hsu* 
aan who describes the ring behavior under density variations. The 
¥ ne erree : . a ‘ . 4 95 » 4 » a ‘ . 
a The line exceeds the surface values by 44-4, 25.0, 8.4, and 2.9 per following brief comments are prompted by personal observations 
—f cent for A = 1, 2, 5, and 10, respectively. The section lift over a period of time. 
curve slope, ao, and the compressibility, 8 = Y%1 — M?, are The direction of growth of the subvortices is found to de 
taken into account by substituting (86Cr,/«) for Crg and (8A /x pend on the relative densities of the ring fluid and the still me 
~| for A, where x = Bao(8)/27 dium and the ring orientation with respect to gravity force —i.e., 
A lifting-surface method is derived which is used to set the up for dye rings less dense than surrounding water, down for 
c limits as 4 + 0 and A — o& for a simple but well-fitted function more dense, and no growth for no density difference. However, 
f \ssume an elliptic load distribution on a straight-load line and the number of such subrings and the speed with which they 
in elliptic chord distribution; and, at the three-quarter chord form is a function of other properties of the ring. Most im 
line, assume that the induced angle of attack is constant span portant is the viscosity of the fluid. Dyes and inks are sus- 
i, | wise and equal to the wing angle of attack. With these condi pensions with non-Newtonian viscosity. If rings are formed by 
tions, a solution of the integral equation of loading is simply a drops of the same volume but different concentrations of ink, 
the initial number of subrings and the number of steps in the 





single solution at the wing plane of symmetry. The results are 


Cra = 27A/[(A/R)E + 2] (1 


cascade process increase with concentration. (The density can 


be adjusted by adding a little alcohol.) Also, using dissolved 
salts as, or with, the marking fluid makes the ring quite active 
Since the instability of the ring seems closely connected to the 
k = (A/4)/W1 + (9A/4)? nonhomogeneous, non-Newtonian fluid, this chain reaction 

may very well be peculiar to this case and not descriptive of the 
As this insta 


mee==| 


where E is the complete elliptic integral of the second kind with 


modulus k, and 


eT 


| 


In the limit conditions as .1 ~ ©, E/k — 1; and, as A > 0, usual transition from laminar to turbulent flow 


E/k — 2/A; thus this Cyq approaches the line-theory value for bility is only manifest when the ring is moving with small trans 
large A and slender-wing theory for small A lational velocity relative to the receiving fluid (not specifically 
For determining a fitted function, assume mentioned in reference 2), the shearing force should be consider 
ably smaller than that occurring in a boundary layer 
The vorticity of the ring fluid probably plays an important 


Liilil 


Cra = 24A/[A + 2f(A)] 


Liit 


- Taking the first and second derivatives of Cyq with respect to A : - 7 . : Mee ; 
a ng . . : ae ee : ‘ Lex x ape hee part in the instability mechanism if the redistribution of vor 
ind likewise of Eq. (1), then, letting A — 0, the derivatives of . : ar : 
, , ; 7 ticity shown by the subvortices is taken as significant If the 
J Eq. (1) require that (a) f(0) = 2, and (b) f’(0) = — (1/2) Also, . ; ; 4 
5 Me , mathematical formulation of such an instability could be made, 
: similarly, for 4 — o, then (c) f(~) = 1, and (d) f’(~) = 0 : ‘ : ma 
- a eh: : it would no doubt add to the knowledge of viscous vortices The 
4 \ssume f(A) = (a; + A)/(a2 + A) Phen condition (a) requires al le ; ; 
1 ‘ Rate: eas similarity of this phenomenon to others such as the sw ellings on 
- that a; = 2a», and (c) is satisfied by the form of the chosen func ie : es aie . 
: a ce ng z ‘ one . the rim of a distended water drop falling in air, the beading of 
. ion. Condition (b) requires that a2 = 2, and condition (d) 1s or — ; ; , A 
| | tisfied } ha eee : TH fA ; ' liquid at the edge of a spinning disk sprayer,' and the breakup 
aso satished Dy the n of the function 1en f(A) = | tT A) . ke . . = . ons 
9 ; 7 — ; Fy ‘ of liquid jets indicate that such vorticity instability effects may 
<+ A) cae ‘ ‘ 
The lift : , share the responsibility with surface tension effects. Recall 
eh ‘urve > g1Ivel \ . ' : re . ‘ - : 
: a Reiner’ credited vorticity with the “teapot effect’’ so long ascribed 


Cre = 247A/{A + 2[(A + 4)/(A + 2)]} to surface tension 
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Fr > 4 stares 4 7 > a ~e > - - - . . - ye ~ es = is i i 
Phe significance of the cascade process to turbulence transition Predictior 
is felt to be due to the interaction of the vortex rings. The EROURED Hg Separatio 
vorticity induced by a tripping device projected into the bound- 2 
ary layer takes time to grow into turbulence, but the details of ‘ EASURED 8 ASSUMING Rgoge® anes M. Rob 
the growth period are hazy. Since organized vorticity is known mS Professor, Det 
to be present at the beginning of the period and the random Rs University ¢ 
distribution of vorticity characteristic of turbulence is the end R; May 6, 1957 
eins ‘ ; ae “qa 2} Ms @ ,REF 2 
result, the possibility of studying viscous vorticity modification 6 on 
and interaction to arrive at a mechanism for transition seems to 8 zo 2 a 
° < =a owe s =z + » 
hold promise 
REFERENCES ee ee ee a a a a flow in 
' 2 3 4 5 6 ? 8 4 0 ration. er 
Sadron, C., Contribution a l'étude experimentale des phenomenes aei “ ; Snenc 
lynamiques aux faibles vitesses, Publ. Sci. et Tech. du Min. de l’Air, No. 22 5 a D. A. spe 
1933 Fic. 1. Ratio of radii of shock wave and body vs. free-streg, | yer may OU 
2 Hsu, C.S., Vortex-Ring Cascade and a Conjectured Model of Transition to Mach Number ut separatio 
Turbulence, Readers’ Forum, Journal of the Aeronautical Sciences, Vol. 24 wuld delay 
No. 4, pp. 318, 314, April, 1957 +e ; 
'’ S 27 
+ Magarvey, R. H., and Taylor, B. W., Free Fall Breakup of Large Drop MEASUREMENTS dition 
J. Appl. Phys., Vol. 27. No. 10, October, 1956 The preset 
‘ Lane, W. R., and Gree ‘ » Mechanics 0 ops ¢ bbles +r . let 
Lane, W R and REED: H.1 The Mechanic f Dre ps and Bu bie lo measure the radius of curvature of the shock wave Rs.t yulent bound 
Surveys in Mechanics, p. 187, Plates 2 and 3; also p. 172: Cambridge Univer eee age z ies regardi 
the Wickes: 1056 position of the shock wave in the vicinity of the axis of symmetry} "0S "8 ' 
Reiner, M., The Teapot Effect. a Problem. Physics Today. Vol. 9. No. 9 was measured on a modified Universal Telereader. The m mndary 1a 
September. 1956 ured points were fitted by both a sphere and parabola. In 
cases the two average values of Rs differed by no more than : \ 
per cent. It is important to take only points ‘sufficiently cloy 
to the vertex of the shock. The criterion for being sufiicient The predic 
> ° ° . ° 
close was that succeeding points should not give consistently j leveloped b 
creasing values of Rs. The results are presented in Fig. 1 wher} existence of 
Rx is the body radius of curvature. sion, Ross. 
The measurements were made on shadowgraphs and interfero.} jye the outer 
-W t i a é ae 
~—— oe ure for Hypersonic grams from three sources: (1) for M = 1.73 and 2.17 shadow} mately relat: 
ri ; 385 i 7 
xisymmetric ow graphs of spherically blunted cones, Rg = 0.385 in. and 0.77)1 proach are ¢ 
in., from the Ballistic Research Laboratory supersonic wind} |ogarithmic | 
R. Sedney tunnel; (2) for M = 3.05, 4.22, and 5.08 shadowgraphs of spheres,} ably depictit 
Ballistic Research Laboratories, Aberdeen Proving Ground, Md Ry = 0.25 in., from the BRL transonic range; (3) for M = 6.66) acted on by s 
9 : eg: . 
April 22, 1957 and 9.4 interferograms of spheres, Rg = 0.25 in., from the BRI For the o 


controlled temperature-pressure range. For the latter the work jaw of Dare 


INTRODUCTION ing medium was Nz at 90°K. For the higher Mach Number | appropriate 


HERE IS GREAT current interest in the hypersonic flow prob interferograms are convenient because the optical distortion j 
7. involving blunt bodies. The flow near the stagnation negligible while the shock wave is still easy'to locate 
point is of particular interest since, for certain kinds of body The shock detachment distance 6 was also measured. Corr ere u is tl 
shapes, the aerodynamic heating is most severe in this region tions for optical distortion were made when necessary utside the 


A number of attempts have been made to represent analytically ) thickness, an 


the axisymmetric flow in this region. Only two of these will be 





DISCUSSION | Ross* * has 
mentioned here. Li and Geiger! have given an approximate red velocit: 
solution and find that a number of their computed results agree In Fig. 1 the measured values of Rs/Rz are given as a function | flat plates I 
with experimental measurements. In particular, the stagnation of M. The four points for the two lowest Mach Numbers a | ew 
point velocity gradient—made suitably nondimensional— is found not joined to the curve because the dissimilarity of the sphericall onde © 
to agree fairly well with experimental results even for Mach blunted cones is noticed at these low Mach Numbers. This i cities 
Numbers as low as 2. However, in arriving at their final results evidenced by the fact that the higher points correspond to t fer predh 
Li and Geiger assume that the radii of curvature of the shock larger value of Rg. The curve, properly extended, should fa ple nN 
wave and body are equal. The purpose of this note is to report below the points for AJ = 1.73 and 2.17 ian wtih 
some measurements which were made of the ratio of these radii Using measured values of 6, the curve obtained by assuming ie 
It is found that if 2.5 7 < 9, where J/ is the free-stream Mach Rs = Rpg + 6 is shown also. It is seen that this approximatior 
Number, this ratio decreases from 1.7 to 1.2 underestimates Rs, and this is noticeable even for the highest 5U9/1 

If the measured ratio is used instead of the assumed value of Mach Numbers ; 
unity, the results of reference 1, say for the stagnation-point Finally, the result from reference 2 is given This analysis here the si 
velocity gradient, will disagree considerably with experimental uses a perturbation method starting from the basic Newton found to a 
results for the lower Mach Numbers. Of course the analysis of flow results. The perturbation parameters are M~? and ( l the factor A 
reference 1 was specifically for hypersonic flows where the ratio of (y+ 1 For the case shown y = 1.4 was used. The agreement | “4ntpressu 
radii is quite close to unity It would appear that, in some man is fairly good for the higher Mach Numbers. It is interesting to | ™etrical floy 
ner, compensating errors were made in the analysis of reference 1 note that agreement over a larger range of M would be obtained | Peunsylvani 
so that agreement at lower Mach Numbers was obtained. Thus if Rs/Rp was assumed to be independent of M/, taking it equalt f Illinois ha 
it is misleading to take this agreement to be verification of the the value for 7 = «a The wall s 
analytical solution 

It is interesting to notice that if one takes the modified New nr ene . 
tonian approximation—i.e., Cp = Cp», cos? 8, where 8 is the where Rg is 
angle between the normal to the body and the free-stream direc Li, T. Y., and Geiger, R. E., Stagnation Point of a Blunt Body in H JaTameter 
tion together with 7 = 1.4—then the Bernoulli equation gives a sonic Flow, Journal of the Aeronautical Sciences, Vol. 24, No. i, pp. 29-% 68D rh 
value of the stagnation-point velocity gradient which differs January, 1957 : veers? roximate re 
ae eal . s y 2 Chester, W., Supersonic Flow Past a Bluff Body with a Deta 
insignificantly from that of reference 1 and thus agrees with Pust'll. Meisemmaitcd Body. Gonsust of Wid Seckenics. Cel 
experimental data even at the lower Mach Numbers p. 490. November, 1956 
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INTRODUCTION 








\E OF THB MAJOR PROBLEMS Of turbulent boundary-layer 
() flow in adverse pressure gradients is the prediction of sep 
sation. Several papers have recently discussed this question 
p. A. Spence? notes that the velocity outside the boundary 
wer may only decrease in a ratio of seven tenths to a half with 
ut separation and that an increase in overall Reynolds Number 
wuld delay separation. H. Schuh* * relates the separation 
ditions to the initial value of the shape parameter H 
The present ce ummunication indicates another method of tur 
lent boundary-layer analysis which leads to qualitative indica 
ions regarding the occurrence of separation in two-dimensional 


undarv layers subjected to adverse pressure gradients 


MsTHOD OF BOUNDARY-LAYER ANALYSIS 


The prediction of separation is based on an analysis method 
leveloped by D Ross 4.6 
xistence of an inner, or wall, and an outer velocity distribution 


Following the observation of the 


region, Ross employed a separate form parameter D to characte1 
Since the flow in the outer region is not inti 
at the wall, the advantages of this ap 


ize the outer flow 
mately related to that 
For the inner region, Ross uses the 


proach are considerable 


garithmic wall-type law which is commonly accepted as suit 





bly depicting the velocity distribution, even when the flow is 
ed on by strong pressure gradients 

For the outer flow region, the three-halves power-deficiency 
iw of Darey (1858) and Prandtl (1925) is found to be most 


:ppropriate 


| (u/U) = Di y/6 )3/? (1 


vhere u is the velocity at wall distance y, U is the velocity just 


utside the boundary layer, 6 is a boundary-layer disturbance 


) thickness, and D is a parameter descriptive of outer region alone 


Ross’ § has found good agreement between this law and meas 
red velocity profiles under a wide variety of conditions. For 
fat plates he finds D S.4°N ‘¢;/2 and for smooth (and rough* 
ips D = 5.1 'C 2 


is thought to be due to intermittancy 


The differences in D for these two cases 
at the outer edge of the 
joundary layer 

For predicting the behavior of turbulent boundary iayers under 
various external conditions, Ross* ® develops an outer flow pa 
rameter which is a function of distance but not local conditions, 


nus, 


function (A 1(.X 


bly Di 


where the subscript 0 refers to the initial conditions. Eq. (2 


sfound to apply to plane two-dimensional flows of all kinds with 
12V cy) computed for the initial con 
One might expect that axisym 
The 


ind recently at the University 


the factor A 0.025 (1 
tant-pressure boundary laver 
>» 


netrical flows would also follow Eq. (2), but studies at 


Peunsylvania State University” 8 
f Illinois have indicated otherwise 


rhe wall shear stress coefficient is given by Ross’ approach as 


+ 5 log (yRo)|? 3) 
where Rg is the momentum thickness Reynolds Number and the 
arameter y is given by the empirical relation y = 0.90 (1 

68D rhe conventionai shape parameter is given by the ap 


proximate relation illustrated in Fig. 1, 


H = 082+ D+403D4 + 4VG/2 (4) 
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Fic. 1. Relation between shape parameters H and D for various 


local skin-friction coefficients 


IN PLANE Two-DIMENSIONAL FLOW Past SMOOTH 


PLATES 


SEPARATION 


With the aid of a simple empirical expression for the rate of 
boundary-layer growth, developed by Ross and Robertson,’ Eq 
(2) permits one to make quantitative predictions for the occur 
rence separation in two-dimensional flows along smooth plates 
Ross finds that the outer shape parameter D must have a value 


of about 1.3 for separation. Eq. (2) is rewritten as 


(6/0)/D = {A[(X Xo)/00)| + [(6/0)/D]o} (U09/ L'e0 5 
where @ is the boundary-laver momentum thickness. Within 
certain limitations, ® 

0/09 = (Ue/U)?*¢ (6 
where G is a unique function of the initial value of Re. One thus 
obtains the relation 

f = (6/0)/D = [A(A&X/0.) + fo] (U/ Ue)? * ” 


where f = (6/0)/D = function (D) 

Ross finds the linear relation of Eq. (2) to hold for AX 
50. At the initial station the conversion from AX /@, to AX /d9 
is a unique function of Reynolds Number, and we get the follow 


50 Ss 


ing relation for the allowable velocity reduction without sepa 
ration: 

fl 
$5/[A "(AX /b0) + fo]; 1" (8 


U;/Uo i fs/[A (AX /b9) 4 
where the subscript s is used to denote conditions at separation 
The value of / 5 corresponds to the separation condition D 
1.3 


initial conditions as represented by 


and fo are functions of the 


the Reynolds Number Rg 


The other parameters cS A. 


rhe solution does depend upon how fast separation is approached 

i.e., on AX /ébg 
verse pressure gradient range 
of Reynolds Number of interest is shown in Fig. 2. The effect 
of high initial Reynolds Numbers and slower rates of velocity 


but apparently not on the details of the ad 


The solution of Eq. (8) for the 


re clearly evident 


deceleration in delaying separation ; 
Comparison of the terminal conditions with those found ex 
Excellent agreement is 


perimentally is also shown in Fig 


Newman” while only 


Schubauer and P. S 


ipparent with the experiments by B. G 
fair agreement appears with that of G. B 


Klebanoff. ! 


EFFECT OF SURFACE ROUGHNESS ON SEPARATION 


The analysis just presented has been for plane two-dimensional 
turbulent Actual 
surfaces are sometimes rough, and one is therefore led to the 
A review 


boundary-layer flow past a smooth surface 


question of whether this helps or inhibits separation 
of our knowledge of the effect of roughness on turbulent boundary 
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Fic. 2. Free-stream velocity reduction ratio for separation 


~ 


layers® '? indicates that the effect of roughness is limited to the 
inncr, or wall, region. Except as the skin friction enters into D, 
the outer flow is independent of surface roughness in the zero 
pressure gradient case. The outer shape parameter D is found 
to be given by the same relation in terms of cy as for smooth 
surfaces. 

Considering cy, rather than Rg, as the basic parameter defining 
the flow, it appears that we can simply transfer our analysis to 
rough surfaces with no further ado. The scale for cy is thus 
added to the Ro abscissa in Fig. 2 to include rough surfaces in the 
prediction. 

The resultant prediction appears reasonable although it seems 
a bit surprising that roughness effects enter in only through the 


initial conditions. 
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IT Eg. (4) y should read 


Boa 1/2 1/2 
ez)y 
n+ 1 { 
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